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