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