QGIS API Documentation  2.99.0-Master (5753576)
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;
123  const char *myTail;
124  sqlite3_stmt *myPreparedStatement;
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::srsDbFilePath().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.isEmpty() )
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.isEmpty() )
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.isEmpty() )
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:148
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.
#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:193
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:295
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
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: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:143
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:163
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.
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.
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:147