QGIS API Documentation  2.99.0-Master (6a61179)
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 <QObject>
21 
22 #include "qgis.h"
23 #include "qgspoint.h"
24 #include "qgscoordinatetransform.h"
26 #include "qgsgeometry.h"
27 #include "qgsgeometrycollection.h"
28 #include "qgsdistancearea.h"
29 #include "qgsapplication.h"
30 #include "qgslogger.h"
31 #include "qgsmessagelog.h"
32 #include "qgsmultisurface.h"
33 #include "qgswkbptr.h"
34 #include "qgslinestring.h"
35 #include "qgspolygon.h"
36 #include "qgssurface.h"
37 #include "qgsunittypes.h"
38 #include "qgscsexception.h"
39 
40 // MSVC compiler doesn't have defined M_PI in math.h
41 #ifndef M_PI
42 #define M_PI 3.14159265358979323846
43 #endif
44 
45 #define DEG2RAD(x) ((x)*M_PI/180)
46 
47 
49 {
50  // init with default settings
51  mEllipsoidalMode = false;
52  setSourceCrs( GEOCRS_ID ); // WGS 84
54 }
55 
56 
59 {
60  _copy( origDA );
61 }
62 
64 {
65 }
66 
69 {
70  if ( this == & origDA )
71  {
72  // Do not copy unto self
73  return *this;
74  }
75  _copy( origDA );
76  return *this;
77 }
78 
80 void QgsDistanceArea::_copy( const QgsDistanceArea & origDA )
81 {
82  mEllipsoidalMode = origDA.mEllipsoidalMode;
83  mEllipsoid = origDA.mEllipsoid;
84  mSemiMajor = origDA.mSemiMajor;
85  mSemiMinor = origDA.mSemiMinor;
86  mInvFlattening = origDA.mInvFlattening;
87  // Some calculations and trig. Should not be TOO time consuming.
88  // Alternatively we could copy the temp vars?
90  mCoordTransform = origDA.mCoordTransform;
91 }
92 
94 {
95  mEllipsoidalMode = flag;
96 }
97 
99 {
100  return mEllipsoidalMode && mEllipsoid != GEO_NONE;
101 }
102 
104 {
106  mCoordTransform.setSourceCrs( srcCRS );
107 }
108 
110 {
111  mCoordTransform.setSourceCrs( srcCRS );
112 }
113 
114 void QgsDistanceArea::setSourceAuthId( const QString& authId )
115 {
117  mCoordTransform.setSourceCrs( srcCRS );
118 }
119 
121 {
122  QString radius, parameter2;
123  //
124  // SQLITE3 stuff - get parameters for selected ellipsoid
125  //
126  sqlite3 *myDatabase;
127  const char *myTail;
128  sqlite3_stmt *myPreparedStatement;
129  int myResult;
130 
131  // Shortcut if ellipsoid is none.
132  if ( ellipsoid == GEO_NONE )
133  {
134  mEllipsoid = GEO_NONE;
135  return true;
136  }
137 
138  // Check if we have a custom projection, and set from text string.
139  // Format is "PARAMETER:<semi-major axis>:<semi minor axis>
140  // Numbers must be with (optional) decimal point and no other separators (C locale)
141  // Distances in meters. Flattening is calculated.
142  if ( ellipsoid.startsWith( QLatin1String( "PARAMETER" ) ) )
143  {
144  QStringList paramList = ellipsoid.split( ':' );
145  bool semiMajorOk, semiMinorOk;
146  double semiMajor = paramList[1].toDouble( & semiMajorOk );
147  double semiMinor = paramList[2].toDouble( & semiMinorOk );
148  if ( semiMajorOk && semiMinorOk )
149  {
150  return setEllipsoid( semiMajor, semiMinor );
151  }
152  else
153  {
154  return false;
155  }
156  }
157 
158  // Continue with PROJ.4 list of ellipsoids.
159 
160  //check the db is available
161  myResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().data(), &myDatabase, SQLITE_OPEN_READONLY, nullptr );
162  if ( myResult )
163  {
164  QgsMessageLog::logMessage( QObject::tr( "Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
165  // XXX This will likely never happen since on open, sqlite creates the
166  // database if it does not exist.
167  return false;
168  }
169  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
170  QString mySql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + '\'';
171  myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
172  // XXX Need to free memory from the error msg if one is set
173  if ( myResult == SQLITE_OK )
174  {
175  if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
176  {
177  radius = QString( reinterpret_cast< const char * >( sqlite3_column_text( myPreparedStatement, 0 ) ) );
178  parameter2 = QString( reinterpret_cast< const char * >( sqlite3_column_text( myPreparedStatement, 1 ) ) );
179  }
180  }
181  // close the sqlite3 statement
182  sqlite3_finalize( myPreparedStatement );
183  sqlite3_close( myDatabase );
184 
185  // row for this ellipsoid wasn't found?
186  if ( radius.isEmpty() || parameter2.isEmpty() )
187  {
188  QgsDebugMsg( QString( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
189  return false;
190  }
191 
192  // get major semiaxis
193  if ( radius.left( 2 ) == QLatin1String( "a=" ) )
194  mSemiMajor = radius.midRef( 2 ).toDouble();
195  else
196  {
197  QgsDebugMsg( QString( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
198  return false;
199  }
200 
201  // get second parameter
202  // one of values 'b' or 'f' is in field parameter2
203  // second one must be computed using formula: invf = a/(a-b)
204  if ( parameter2.left( 2 ) == QLatin1String( "b=" ) )
205  {
206  mSemiMinor = parameter2.midRef( 2 ).toDouble();
207  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
208  }
209  else if ( parameter2.left( 3 ) == QLatin1String( "rf=" ) )
210  {
211  mInvFlattening = parameter2.midRef( 3 ).toDouble();
212  mSemiMinor = mSemiMajor - ( mSemiMajor / mInvFlattening );
213  }
214  else
215  {
216  QgsDebugMsg( QString( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
217  return false;
218  }
219 
220  QgsDebugMsg( QString( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
221 
222 
223  // get spatial ref system for ellipsoid
224  QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
226  //TODO: createFromProj4 used to save to the user database any new CRS
227  // this behavior was changed in order to separate creation and saving.
228  // Not sure if it necessary to save it here, should be checked by someone
229  // familiar with the code (should also give a more descriptive name to the generated CRS)
230  if ( destCRS.srsid() == 0 )
231  {
232  QString myName = QStringLiteral( " * %1 (%2)" )
233  .arg( QObject::tr( "Generated CRS", "A CRS automatically generated from layer info get this prefix for description" ),
234  destCRS.toProj4() );
235  destCRS.saveAsUserCrs( myName );
236  }
237  //
238 
239  // set transformation from project CRS to ellipsoid coordinates
240  mCoordTransform.setDestinationCrs( destCRS );
241 
242  mEllipsoid = ellipsoid;
243 
244  // precalculate some values for area calculations
245  computeAreaInit();
246 
247  return true;
248 }
249 
251 // Inverse flattening is calculated with invf = a/(a-b)
252 // Also, b = a-(a/invf)
253 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
254 {
255  mEllipsoid = QStringLiteral( "PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
256  mSemiMajor = semiMajor;
257  mSemiMinor = semiMinor;
258  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
259 
260  computeAreaInit();
261 
262  return true;
263 }
264 
265 double QgsDistanceArea::measure( const QgsAbstractGeometry* geomV2, MeasureType type ) const
266 {
267  if ( !geomV2 )
268  {
269  return 0.0;
270  }
271 
272  int geomDimension = geomV2->dimension();
273  if ( geomDimension <= 0 )
274  {
275  return 0.0;
276  }
277 
278  MeasureType measureType = type;
279  if ( measureType == Default )
280  {
281  measureType = ( geomDimension == 1 ? Length : Area );
282  }
283 
284  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
285  {
286  //no transform required
287  if ( measureType == Length )
288  {
289  return geomV2->length();
290  }
291  else
292  {
293  return geomV2->area();
294  }
295  }
296  else
297  {
298  //multigeom is sum of measured parts
299  const QgsGeometryCollection* collection = dynamic_cast<const QgsGeometryCollection*>( geomV2 );
300  if ( collection )
301  {
302  double sum = 0;
303  for ( int i = 0; i < collection->numGeometries(); ++i )
304  {
305  sum += measure( collection->geometryN( i ), measureType );
306  }
307  return sum;
308  }
309 
310  if ( measureType == Length )
311  {
312  const QgsCurve* curve = dynamic_cast<const QgsCurve*>( geomV2 );
313  if ( !curve )
314  {
315  return 0.0;
316  }
317 
318  QgsLineString* lineString = curve->curveToLine();
319  double length = measureLine( lineString );
320  delete lineString;
321  return length;
322  }
323  else
324  {
325  const QgsSurface* surface = dynamic_cast<const QgsSurface*>( geomV2 );
326  if ( !surface )
327  return 0.0;
328 
329  QgsPolygonV2* polygon = surface->surfaceToPolygon();
330 
331  double area = 0;
332  const QgsCurve* outerRing = polygon->exteriorRing();
333  area += measurePolygon( outerRing );
334 
335  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
336  {
337  const QgsCurve* innerRing = polygon->interiorRing( i );
338  area -= measurePolygon( innerRing );
339  }
340  delete polygon;
341  return area;
342  }
343  }
344 }
345 
346 double QgsDistanceArea::measureArea( const QgsGeometry* geometry ) const
347 {
348  if ( !geometry )
349  return 0.0;
350 
351  const QgsAbstractGeometry* geomV2 = geometry->geometry();
352  return measure( geomV2, Area );
353 }
354 
355 double QgsDistanceArea::measureArea( const QgsGeometry& geometry ) const
356 {
357  if ( geometry.isEmpty() )
358  return 0.0;
359 
360  const QgsAbstractGeometry* geomV2 = geometry.geometry();
361  return measure( geomV2, Area );
362 }
363 
364 double QgsDistanceArea::measureLength( const QgsGeometry* geometry ) const
365 {
366  if ( !geometry )
367  return 0.0;
368 
369  const QgsAbstractGeometry* geomV2 = geometry->geometry();
370  return measure( geomV2, Length );
371 }
372 
373 double QgsDistanceArea::measureLength( const QgsGeometry& geometry ) const
374 {
375  if ( geometry.isEmpty() )
376  return 0.0;
377 
378  const QgsAbstractGeometry* geomV2 = geometry.geometry();
379  return measure( geomV2, Length );
380 }
381 
382 double QgsDistanceArea::measurePerimeter( const QgsGeometry* geometry ) const
383 {
384  if ( !geometry )
385  return 0.0;
386 
387  return measurePerimeter( *geometry );
388 }
389 
390 double QgsDistanceArea::measurePerimeter( const QgsGeometry& geometry ) const
391 {
392  if ( geometry.isEmpty() )
393  return 0.0;
394 
395  const QgsAbstractGeometry* geomV2 = geometry.geometry();
396  if ( !geomV2 || geomV2->dimension() < 2 )
397  {
398  return 0.0;
399  }
400 
401  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
402  {
403  return geomV2->perimeter();
404  }
405 
406  //create list with (single) surfaces
407  QList< const QgsSurface* > surfaces;
408  const QgsSurface* surf = dynamic_cast<const QgsSurface*>( geomV2 );
409  if ( surf )
410  {
411  surfaces.append( surf );
412  }
413  const QgsMultiSurface* multiSurf = dynamic_cast<const QgsMultiSurface*>( geomV2 );
414  if ( multiSurf )
415  {
416  surfaces.reserve(( surf ? 1 : 0 ) + multiSurf->numGeometries() );
417  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
418  {
419  surfaces.append( static_cast<const QgsSurface*>( multiSurf->geometryN( i ) ) );
420  }
421  }
422 
423  double length = 0;
424  QList<const QgsSurface*>::const_iterator surfaceIt = surfaces.constBegin();
425  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
426  {
427  if ( !*surfaceIt )
428  {
429  continue;
430  }
431 
432  QgsPolygonV2* poly = ( *surfaceIt )->surfaceToPolygon();
433  const QgsCurve* outerRing = poly->exteriorRing();
434  if ( outerRing )
435  {
436  length += measure( outerRing );
437  }
438  int nInnerRings = poly->numInteriorRings();
439  for ( int i = 0; i < nInnerRings; ++i )
440  {
441  length += measure( poly->interiorRing( i ) );
442  }
443  delete poly;
444  }
445  return length;
446 }
447 
448 double QgsDistanceArea::measureLine( const QgsCurve* curve ) const
449 {
450  if ( !curve )
451  {
452  return 0.0;
453  }
454 
455  QgsPointSequence linePointsV2;
456  QList<QgsPoint> linePoints;
457  curve->points( linePointsV2 );
458  QgsGeometry::convertPointList( linePointsV2, linePoints );
459  return measureLine( linePoints );
460 }
461 
462 double QgsDistanceArea::measureLine( const QList<QgsPoint> &points ) const
463 {
464  if ( points.size() < 2 )
465  return 0;
466 
467  double total = 0;
468  QgsPoint p1, p2;
469 
470  try
471  {
472  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
473  p1 = mCoordTransform.transform( points[0] );
474  else
475  p1 = points[0];
476 
477  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
478  {
479  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
480  {
481  p2 = mCoordTransform.transform( *i );
482  total += computeDistanceBearing( p1, p2 );
483  }
484  else
485  {
486  p2 = *i;
487  total += measureLine( p1, p2 );
488  }
489 
490  p1 = p2;
491  }
492 
493  return total;
494  }
495  catch ( QgsCsException &cse )
496  {
497  Q_UNUSED( cse );
498  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
499  return 0.0;
500  }
501 
502 }
503 
504 double QgsDistanceArea::measureLine( const QgsPoint &p1, const QgsPoint &p2 ) const
505 {
507  return measureLine( p1, p2, units );
508 }
509 
510 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2, QgsUnitTypes::DistanceUnit& units ) const
511 {
512  double result;
513  units = mCoordTransform.sourceCrs().mapUnits();
514 
515  try
516  {
517  QgsPoint pp1 = p1, pp2 = p2;
518 
519  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
520  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
521  {
523  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
524  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform.sourceCrs().toProj4() ), 4 );
525  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj4() ), 4 );
526  pp1 = mCoordTransform.transform( p1 );
527  pp2 = mCoordTransform.transform( p2 );
528  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
529  result = computeDistanceBearing( pp1, pp2 );
530  }
531  else
532  {
533  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
534  result = computeDistanceFlat( p1, p2 );
535  }
536  }
537  catch ( QgsCsException &cse )
538  {
539  Q_UNUSED( cse );
540  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
541  result = 0.0;
542  }
543  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
544  return result;
545 }
546 
548 {
549  return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();
550 }
551 
553 {
555  QgsUnitTypes::distanceToAreaUnit( mCoordTransform.sourceCrs().mapUnits() );
556 }
557 
558 QgsConstWkbPtr QgsDistanceArea::measurePolygon( QgsConstWkbPtr wkbPtr, double* area, double* perimeter, bool hasZptr ) const
559 {
560  if ( !wkbPtr )
561  {
562  QgsDebugMsg( "no feature to measure" );
563  return wkbPtr;
564  }
565 
566  wkbPtr.readHeader();
567 
568  // get number of rings in the polygon
569  int numRings;
570  wkbPtr >> numRings;
571 
572  if ( numRings == 0 )
573  {
574  QgsDebugMsg( "no rings to measure" );
575  return QgsConstWkbPtr( nullptr, 0 );
576  }
577 
578  // Set pointer to the first ring
579  QList<QgsPoint> points;
580  QgsPoint pnt;
581  double x, y;
582  if ( area )
583  *area = 0;
584  if ( perimeter )
585  *perimeter = 0;
586 
587  try
588  {
589  for ( int idx = 0; idx < numRings; idx++ )
590  {
591  int nPoints;
592  wkbPtr >> nPoints;
593 
594  // Extract the points from the WKB and store in a pair of
595  // vectors.
596  for ( int jdx = 0; jdx < nPoints; jdx++ )
597  {
598  wkbPtr >> x >> y;
599  if ( hasZptr )
600  {
601  // totally ignore Z value
602  wkbPtr += sizeof( double );
603  }
604 
605  pnt = QgsPoint( x, y );
606 
607  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
608  {
609  pnt = mCoordTransform.transform( pnt );
610  }
611  points.append( pnt );
612  }
613 
614  if ( points.size() > 2 )
615  {
616  if ( area )
617  {
618  double areaTmp = computePolygonArea( points );
619  if ( idx == 0 )
620  {
621  // exterior ring
622  *area += areaTmp;
623  }
624  else
625  {
626  *area -= areaTmp; // interior rings
627  }
628  }
629 
630  if ( perimeter )
631  {
632  if ( idx == 0 )
633  {
634  // exterior ring
635  *perimeter += computeDistance( points );
636  }
637  }
638  }
639 
640  points.clear();
641  }
642  }
643  catch ( QgsCsException &cse )
644  {
645  Q_UNUSED( cse );
646  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
647  }
648 
649  return wkbPtr;
650 }
651 
652 double QgsDistanceArea::measurePolygon( const QgsCurve* curve ) const
653 {
654  if ( !curve )
655  {
656  return 0.0;
657  }
658 
659  QgsPointSequence linePointsV2;
660  curve->points( linePointsV2 );
661  QList<QgsPoint> linePoints;
662  QgsGeometry::convertPointList( linePointsV2, linePoints );
663  return measurePolygon( linePoints );
664 }
665 
666 
667 double QgsDistanceArea::measurePolygon( const QList<QgsPoint>& points ) const
668 {
669  try
670  {
671  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
672  {
673  QList<QgsPoint> pts;
674  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
675  {
676  pts.append( mCoordTransform.transform( *i ) );
677  }
678  return computePolygonArea( pts );
679  }
680  else
681  {
682  return computePolygonArea( points );
683  }
684  }
685  catch ( QgsCsException &cse )
686  {
687  Q_UNUSED( cse );
688  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
689  return 0.0;
690  }
691 }
692 
693 
694 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 ) const
695 {
696  QgsPoint pp1 = p1, pp2 = p2;
697  double bearing;
698 
699  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
700  {
701  pp1 = mCoordTransform.transform( p1 );
702  pp2 = mCoordTransform.transform( p2 );
703  computeDistanceBearing( pp1, pp2, &bearing );
704  }
705  else //compute simple planar azimuth
706  {
707  double dx = p2.x() - p1.x();
708  double dy = p2.y() - p1.y();
709  bearing = atan2( dx, dy );
710  }
711 
712  return bearing;
713 }
714 
715 
717 // distance calculation
718 
720  const QgsPoint& p1, const QgsPoint& p2,
721  double* course1, double* course2 ) const
722 {
723  if ( qgsDoubleNear( p1.x(), p2.x() ) && qgsDoubleNear( p1.y(), p2.y() ) )
724  return 0;
725 
726  // ellipsoid
727  double a = mSemiMajor;
728  double b = mSemiMinor;
729  double f = 1 / mInvFlattening;
730 
731  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
732  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
733 
734  double L = p2_lon - p1_lon;
735  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
736  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
737  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
738  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
739  double lambda = L;
740  double lambdaP = 2 * M_PI;
741 
742  double sinLambda = 0;
743  double cosLambda = 0;
744  double sinSigma = 0;
745  double cosSigma = 0;
746  double sigma = 0;
747  double alpha = 0;
748  double cosSqAlpha = 0;
749  double cos2SigmaM = 0;
750  double C = 0;
751  double tu1 = 0;
752  double tu2 = 0;
753 
754  int iterLimit = 20;
755  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
756  {
757  sinLambda = sin( lambda );
758  cosLambda = cos( lambda );
759  tu1 = ( cosU2 * sinLambda );
760  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
761  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
762  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
763  sigma = atan2( sinSigma, cosSigma );
764  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
765  cosSqAlpha = cos( alpha ) * cos( alpha );
766  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
767  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
768  lambdaP = lambda;
769  lambda = L + ( 1 - C ) * f * sin( alpha ) *
770  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
771  }
772 
773  if ( iterLimit == 0 )
774  return -1; // formula failed to converge
775 
776  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
777  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
778  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
779  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
780  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
781  double s = b * A * ( sigma - deltaSigma );
782 
783  if ( course1 )
784  {
785  *course1 = atan2( tu1, tu2 );
786  }
787  if ( course2 )
788  {
789  // PI is added to return azimuth from P2 to P1
790  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
791  }
792 
793  return s;
794 }
795 
796 double QgsDistanceArea::computeDistanceFlat( const QgsPoint& p1, const QgsPoint& p2 ) const
797 {
798  return sqrt(( p2.x() - p1.x() ) * ( p2.x() - p1.x() ) + ( p2.y() - p1.y() ) * ( p2.y() - p1.y() ) );
799 }
800 
801 double QgsDistanceArea::computeDistance( const QList<QgsPoint>& points ) const
802 {
803  if ( points.size() < 2 )
804  return 0;
805 
806  double total = 0;
807  QgsPoint p1, p2;
808 
809  try
810  {
811  p1 = points[0];
812 
813  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
814  {
815  p2 = *i;
816  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
817  {
818  total += computeDistanceBearing( p1, p2 );
819  }
820  else
821  {
822  total += computeDistanceFlat( p1, p2 );
823  }
824 
825  p1 = p2;
826  }
827 
828  return total;
829  }
830  catch ( QgsCsException &cse )
831  {
832  Q_UNUSED( cse );
833  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
834  return 0.0;
835  }
836 }
837 
838 
839 
841 // stuff for measuring areas - copied from GRASS
842 // don't know how does it work, but it's working .)
843 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
844 
845 double QgsDistanceArea::getQ( double x ) const
846 {
847  double sinx, sinx2;
848 
849  sinx = sin( x );
850  sinx2 = sinx * sinx;
851 
852  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
853 }
854 
855 
856 double QgsDistanceArea::getQbar( double x ) const
857 {
858  double cosx, cosx2;
859 
860  cosx = cos( x );
861  cosx2 = cosx * cosx;
862 
863  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
864 }
865 
866 
868 {
869  //don't try to perform calculations if no ellipsoid
870  if ( mEllipsoid == GEO_NONE )
871  {
872  return;
873  }
874 
875  double a2 = ( mSemiMajor * mSemiMajor );
876  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
877  double e4, e6;
878 
879  m_TwoPI = M_PI + M_PI;
880 
881  e4 = e2 * e2;
882  e6 = e4 * e2;
883 
884  m_AE = a2 * ( 1 - e2 );
885 
886  m_QA = ( 2.0 / 3.0 ) * e2;
887  m_QB = ( 3.0 / 5.0 ) * e4;
888  m_QC = ( 4.0 / 7.0 ) * e6;
889 
890  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
891  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
892  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
893  m_QbarD = ( 4.0 / 49.0 ) * e6;
894 
895  m_Qp = getQ( M_PI / 2 );
896  m_E = 4 * M_PI * m_Qp * m_AE;
897  if ( m_E < 0.0 )
898  m_E = -m_E;
899 }
900 
901 
902 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint>& points ) const
903 {
904  if ( points.isEmpty() )
905  {
906  return 0;
907  }
908 
909  double x1, y1, x2, y2, dx, dy;
910  double Qbar1, Qbar2;
911  double area;
912 
913  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
914  if (( ! mEllipsoidalMode ) || ( mEllipsoid == GEO_NONE ) )
915  {
916  return computePolygonFlatArea( points );
917  }
918  int n = points.size();
919  x2 = DEG2RAD( points[n-1].x() );
920  y2 = DEG2RAD( points[n-1].y() );
921  Qbar2 = getQbar( y2 );
922 
923  area = 0.0;
924 
925  for ( int i = 0; i < n; i++ )
926  {
927  x1 = x2;
928  y1 = y2;
929  Qbar1 = Qbar2;
930 
931  x2 = DEG2RAD( points[i].x() );
932  y2 = DEG2RAD( points[i].y() );
933  Qbar2 = getQbar( y2 );
934 
935  if ( x1 > x2 )
936  while ( x1 - x2 > M_PI )
937  x2 += m_TwoPI;
938  else if ( x2 > x1 )
939  while ( x2 - x1 > M_PI )
940  x1 += m_TwoPI;
941 
942  dx = x2 - x1;
943  area += dx * ( m_Qp - getQ( y2 ) );
944 
945  dy = y2 - y1;
946  if ( !qgsDoubleNear( dy, 0.0 ) )
947  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
948  }
949  if (( area *= m_AE ) < 0.0 )
950  area = -area;
951 
952  /* kludge - if polygon circles the south pole the area will be
953  * computed as if it cirlced the north pole. The correction is
954  * the difference between total surface area of the earth and
955  * the "north pole" area.
956  */
957  if ( area > m_E )
958  area = m_E;
959  if ( area > m_E / 2 )
960  area = m_E - area;
961 
962  return area;
963 }
964 
965 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint>& points ) const
966 {
967  // Normal plane area calculations.
968  double area = 0.0;
969  int i, size;
970 
971  size = points.size();
972 
973  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
974  for ( i = 0; i < size; i++ )
975  {
976  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
977  // Using '% size', so that we always end with the starting point
978  // and thus close the polygon.
979  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
980  }
981  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
982  area = area / 2.0;
983  return qAbs( area ); // All areas are positive!
984 }
985 
986 QString QgsDistanceArea::formatDistance( double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit )
987 {
988  QString unitLabel;
989 
990  switch ( unit )
991  {
993  if ( keepBaseUnit || qAbs( distance ) == 0.0 )
994  {
995  unitLabel = QObject::tr( " m" );
996  }
997  else if ( qAbs( distance ) > 1000.0 )
998  {
999  unitLabel = QObject::tr( " km" );
1000  distance = distance / 1000;
1001  }
1002  else if ( qAbs( distance ) < 0.01 )
1003  {
1004  unitLabel = QObject::tr( " mm" );
1005  distance = distance * 1000;
1006  }
1007  else if ( qAbs( distance ) < 0.1 )
1008  {
1009  unitLabel = QObject::tr( " cm" );
1010  distance = distance * 100;
1011  }
1012  else
1013  {
1014  unitLabel = QObject::tr( " m" );
1015  }
1016  break;
1017 
1019  if ( keepBaseUnit || qAbs( distance ) >= 1.0 )
1020  {
1021  unitLabel = QObject::tr( " km" );
1022  }
1023  else
1024  {
1025  unitLabel = QObject::tr( " m" );
1026  distance = distance * 1000;
1027  }
1028  break;
1029 
1031  if ( qAbs( distance ) <= 5280.0 || keepBaseUnit )
1032  {
1033  unitLabel = QObject::tr( " ft" );
1034  }
1035  else
1036  {
1037  unitLabel = QObject::tr( " mi" );
1038  distance /= 5280.0;
1039  }
1040  break;
1041 
1043  if ( qAbs( distance ) <= 1760.0 || keepBaseUnit )
1044  {
1045  unitLabel = QObject::tr( " yd" );
1046  }
1047  else
1048  {
1049  unitLabel = QObject::tr( " mi" );
1050  distance /= 1760.0;
1051  }
1052  break;
1053 
1055  if ( qAbs( distance ) >= 1.0 || keepBaseUnit )
1056  {
1057  unitLabel = QObject::tr( " mi" );
1058  }
1059  else
1060  {
1061  unitLabel = QObject::tr( " ft" );
1062  distance *= 5280.0;
1063  }
1064  break;
1065 
1067  unitLabel = QObject::tr( " NM" );
1068  break;
1069 
1071 
1072  if ( qAbs( distance ) == 1.0 )
1073  unitLabel = QObject::tr( " degree" );
1074  else
1075  unitLabel = QObject::tr( " degrees" );
1076  break;
1077 
1079  unitLabel.clear();
1080  break;
1081  default:
1082  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( unit ) );
1083  break;
1084  }
1085 
1086  return QStringLiteral( "%L1%2" ).arg( distance, 0, 'f', decimals ).arg( unitLabel );
1087 }
1088 
1089 QString QgsDistanceArea::formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit )
1090 {
1091  QString unitLabel;
1092 
1093  switch ( unit )
1094  {
1096  {
1097  if ( keepBaseUnit )
1098  {
1099  unitLabel = QObject::trUtf8( " m²" );
1100  }
1102  {
1103  unitLabel = QObject::trUtf8( " km²" );
1105  }
1107  {
1108  unitLabel = QObject::tr( " ha" );
1110  }
1111  else
1112  {
1113  unitLabel = QObject::trUtf8( " m²" );
1114  }
1115  break;
1116  }
1117 
1119  {
1120  unitLabel = QObject::trUtf8( " km²" );
1121  break;
1122  }
1123 
1125  {
1126  if ( keepBaseUnit )
1127  {
1128  unitLabel = QObject::trUtf8( " ft²" );
1129  }
1131  {
1132  unitLabel = QObject::trUtf8( " mi²" );
1134  }
1135  else
1136  {
1137  unitLabel = QObject::trUtf8( " ft²" );
1138  }
1139  break;
1140  }
1141 
1143  {
1144  if ( keepBaseUnit )
1145  {
1146  unitLabel = QObject::trUtf8( " yd²" );
1147  }
1149  {
1150  unitLabel = QObject::trUtf8( " mi²" );
1152  }
1153  else
1154  {
1155  unitLabel = QObject::trUtf8( " yd²" );
1156  }
1157  break;
1158  }
1159 
1161  {
1162  unitLabel = QObject::trUtf8( " mi²" );
1163  break;
1164  }
1165 
1167  {
1168  if ( keepBaseUnit )
1169  {
1170  unitLabel = QObject::trUtf8( " ha" );
1171  }
1173  {
1174  unitLabel = QObject::trUtf8( " km²" );
1176  }
1177  else
1178  {
1179  unitLabel = QObject::trUtf8( " ha" );
1180  }
1181  break;
1182  }
1183 
1185  {
1186  if ( keepBaseUnit )
1187  {
1188  unitLabel = QObject::trUtf8( " ac" );
1189  }
1191  {
1192  unitLabel = QObject::trUtf8( " mi²" );
1194  }
1195  else
1196  {
1197  unitLabel = QObject::trUtf8( " ac" );
1198  }
1199  break;
1200  }
1201 
1203  {
1204  unitLabel = QObject::trUtf8( " nm²" );
1205  break;
1206  }
1207 
1209  {
1210  unitLabel = QObject::tr( " sq.deg." );
1211  break;
1212  }
1213 
1215  {
1216  unitLabel.clear();
1217  break;
1218  }
1219  }
1220 
1221  return QStringLiteral( "%L1%2" ).arg( area, 0, 'f', decimals ).arg( unitLabel );
1222 }
1223 
1225 {
1226  // get the conversion factor between the specified units
1227  QgsUnitTypes::DistanceUnit measureUnits = lengthUnits();
1228  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1229 
1230  double result = length * factorUnits;
1231  QgsDebugMsg( QString( "Converted length of %1 %2 to %3 %4" ).arg( length )
1232  .arg( QgsUnitTypes::toString( measureUnits ) )
1233  .arg( result )
1234  .arg( QgsUnitTypes::toString( toUnits ) ) );
1235  return result;
1236 }
1237 
1239 {
1240  // get the conversion factor between the specified units
1241  QgsUnitTypes::AreaUnit measureUnits = areaUnits();
1242  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1243 
1244  double result = area * factorUnits;
1245  QgsDebugMsg( QString( "Converted area of %1 %2 to %3 %4" ).arg( area )
1246  .arg( QgsUnitTypes::toString( measureUnits ) )
1247  .arg( result )
1248  .arg( QgsUnitTypes::toString( toUnits ) ) );
1249  return result;
1250 }
static void convertPointList(const QList< QgsPoint > &input, QgsPointSequence &output)
Upgrades a point list from QgsPoint to QgsPointV2.
const QgsCurve * interiorRing(int i) const
virtual double length() const
Returns the length of the geometry.
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
double y
Definition: qgspoint.h:116
static QgsCoordinateReferenceSystem fromProj4(const QString &proj4)
Creates a CRS from a proj4 style formatted string.
static Q_INVOKABLE QString toString(QgsUnitTypes::DistanceUnit unit)
Returns a translated string representing a distance unit.
virtual QgsPolygonV2 * surfaceToPolygon() const =0
Get a polygon representation of this surface.
bool saveAsUserCrs(const QString &name)
Save the proj4-string as a custom CRS.
~QgsDistanceArea()
Destructor.
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
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)
QString toProj4() const
Returns a Proj4 string representation of this CRS.
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
Unknown areal unit.
Definition: qgsunittypes.h:75
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:78
double convertLengthMeasurement(double length, QgsUnitTypes::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
bool setEllipsoid(const QString &ellipsoid)
Sets ellipsoid by its acronym.
double convertAreaMeasurement(double area, QgsUnitTypes::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
virtual double area() const
Returns the area of the geometry.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:196
Multi surface geometry collection.
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:75
bool willUseEllipsoid() const
Returns whether calculations will use the ellipsoid.
Polygon geometry type.
Definition: qgspolygon.h:29
int numInteriorRings() const
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from...
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
QgsPoint transform(const QgsPoint &point, TransformDirection direction=ForwardTransform) const
Transform the point from the source CRS to the destination CRS.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
Geometry collection.
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
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:298
static void logMessage(const QString &message, const QString &tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
Degrees, for planar geographic CRS distance measurements.
Definition: qgsunittypes.h:50
bool isEmpty() const
Returns true if the geometry is empty (ie, contains no underlying geometry accessible via geometry)...
double measureLine(const QList< QgsPoint > &points) const
Measures the length of a line with multiple segments.
double computePolygonFlatArea(const QList< QgsPoint > &points) const
double measurePolygon(const QList< QgsPoint > &points) const
measures polygon area
Abstract base class for curved geometry type.
Definition: qgscurve.h:32
Abstract base class for all geometries.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
const QgsCurve * exteriorRing() const
Square degrees, for planar geographic CRS area measurements.
Definition: qgsunittypes.h:74
A class to represent a point.
Definition: qgspoint.h:111
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:129
double computePolygonArea(const QList< QgsPoint > &points) const
calculates area of polygon on ellipsoid algorithm has been taken from GRASS: gis/area_poly1.c
int numGeometries() const
Returns the number of geometries within the collection.
struct sqlite3 sqlite3
double measureArea(const QgsGeometry *geometry) const
Measures the area of a geometry.
QgsDistanceArea & operator=(const QgsDistanceArea &origDA)
Assignment operator.
DistanceUnit
Units of distance.
Definition: qgsunittypes.h:42
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2) const
uses flat / planimetric / Euclidean distance
Unknown distance unit.
Definition: qgsunittypes.h:51
General purpose distance and area calculator.
double measurePerimeter(const QgsGeometry *geometry) const
Measures the perimeter of a polygon geometry.
QgsUnitTypes::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
QString ellipsoid() const
Returns ellipsoid&#39;s acronym.
static QgsCoordinateReferenceSystem fromOgcWmsCrs(const QString &ogcCrs)
Creates a CRS from a given OGC WMS-format Coordinate Reference System string.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
double computeDistanceBearing(const QgsPoint &p1, const QgsPoint &p2, double *course1=nullptr, double *course2=nullptr) const
calculates distance from two points on ellipsoid based on inverse Vincenty&#39;s formulae ...
QgsPolygonV2 * surfaceToPolygon() const override
Get a polygon representation of this surface.
Definition: qgspolygon.cpp:298
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:35
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
void setSourceAuthId(const QString &authid)
sets source spatial reference system by authid
This class represents a coordinate reference system (CRS).
static QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
double measureLength(const QgsGeometry *geometry) const
Measures the length of a geometry.
static QString srsDbFilePath()
Returns the path to the srs.db file.
Custom exception class for Coordinate Reference System related exceptions.
long srsid() const
Returns the internal CRS ID, if available.
Terrestial miles.
Definition: qgsunittypes.h:49
static Q_INVOKABLE double fromUnitToUnitFactor(QgsUnitTypes::DistanceUnit fromUnit, QgsUnitTypes::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
QList< QgsPointV2 > QgsPointSequence
QgsWkbTypes::Type readHeader() const
Definition: qgswkbptr.cpp:53
QgsDistanceArea()
Constructor.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
static Q_INVOKABLE AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, eg meters to square meters.
QgsAbstractGeometry * geometry() const
Returns the underlying geometry store.
AreaUnit
Units of area.
Definition: qgsunittypes.h:64
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
void setEllipsoidalMode(bool flag)
Sets whether coordinates must be projected to ellipsoid before measuring.
double computeDistance(const QList< QgsPoint > &points) const
calculate distance with given coordinates (does not do a transform anymore)
static QgsCoordinateReferenceSystem fromSrsId(long srsId)
Creates a CRS from a specified QGIS SRS ID.
static QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
virtual double perimeter() const
Returns the perimeter of the geometry.
double x
Definition: qgspoint.h:115