QGIS API Documentation  2.15.0-Master (02a0ebe)
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 #include "qgscrscache.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  mCoordTransform = new QgsCoordinateTransform;
53  setSourceCrs( GEOCRS_ID ); // WGS 84
55 }
56 
57 
60  : mCoordTransform( nullptr )
61 {
62  _copy( origDA );
63 }
64 
66 {
67  delete mCoordTransform;
68 }
69 
72 {
73  if ( this == & origDA )
74  {
75  // Do not copy unto self
76  return *this;
77  }
78  _copy( origDA );
79  return *this;
80 }
81 
83 void QgsDistanceArea::_copy( const QgsDistanceArea & origDA )
84 {
85  mEllipsoidalMode = origDA.mEllipsoidalMode;
86  mEllipsoid = origDA.mEllipsoid;
87  mSemiMajor = origDA.mSemiMajor;
88  mSemiMinor = origDA.mSemiMinor;
89  mInvFlattening = origDA.mInvFlattening;
90  // Some calculations and trig. Should not be TOO time consuming.
91  // Alternatively we could copy the temp vars?
93  delete mCoordTransform;
94  mCoordTransform = new QgsCoordinateTransform( origDA.mCoordTransform->sourceCrs(), origDA.mCoordTransform->destCRS() );
95 }
96 
98 {
99  mEllipsoidalMode = flag;
100 }
101 
103 {
104  return mEllipsoidalMode && mEllipsoid != GEO_NONE;
105 }
106 
108 {
110  srcCRS.createFromSrsId( srsid );
111  mCoordTransform->setSourceCrs( srcCRS );
112 }
113 
115 {
116  mCoordTransform->setSourceCrs( srcCRS );
117 }
118 
120 {
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  //TODO: createFromProj4 used to save to the user database any new CRS
232  // this behavior was changed in order to separate creation and saving.
233  // Not sure if it necessary to save it here, should be checked by someone
234  // familiar with the code (should also give a more descriptive name to the generated CRS)
235  if ( destCRS.srsid() == 0 )
236  {
237  QString myName = QString( " * %1 (%2)" )
238  .arg( QObject::tr( "Generated CRS", "A CRS automatically generated from layer info get this prefix for description" ),
239  destCRS.toProj4() );
240  destCRS.saveAsUserCRS( myName );
241  }
242  //
243 
244  // set transformation from project CRS to ellipsoid coordinates
245  mCoordTransform->setDestCRS( destCRS );
246 
247  mEllipsoid = ellipsoid;
248 
249  // precalculate some values for area calculations
250  computeAreaInit();
251 
252  return true;
253 }
254 
256 // Inverse flattening is calculated with invf = a/(a-b)
257 // Also, b = a-(a/invf)
258 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
259 {
260  mEllipsoid = QString( "PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
261  mSemiMajor = semiMajor;
262  mSemiMinor = semiMinor;
263  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
264 
265  computeAreaInit();
266 
267  return true;
268 }
269 
270 double QgsDistanceArea::measure( const QgsAbstractGeometryV2* geomV2, MeasureType type ) const
271 {
272  if ( !geomV2 )
273  {
274  return 0.0;
275  }
276 
277  int geomDimension = geomV2->dimension();
278  if ( geomDimension <= 0 )
279  {
280  return 0.0;
281  }
282 
283  MeasureType measureType = type;
284  if ( measureType == Default )
285  {
286  measureType = ( geomDimension == 1 ? Length : Area );
287  }
288 
289  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
290  {
291  //no transform required
292  if ( measureType == Length )
293  {
294  return geomV2->length();
295  }
296  else
297  {
298  return geomV2->area();
299  }
300  }
301  else
302  {
303  //multigeom is sum of measured parts
304  const QgsGeometryCollectionV2* collection = dynamic_cast<const QgsGeometryCollectionV2*>( geomV2 );
305  if ( collection )
306  {
307  double sum = 0;
308  for ( int i = 0; i < collection->numGeometries(); ++i )
309  {
310  sum += measure( collection->geometryN( i ), measureType );
311  }
312  return sum;
313  }
314 
315  if ( measureType == Length )
316  {
317  const QgsCurveV2* curve = dynamic_cast<const QgsCurveV2*>( geomV2 );
318  if ( !curve )
319  {
320  return 0.0;
321  }
322 
323  QgsLineStringV2* lineString = curve->curveToLine();
324  double length = measureLine( lineString );
325  delete lineString;
326  return length;
327  }
328  else
329  {
330  const QgsSurfaceV2* surface = dynamic_cast<const QgsSurfaceV2*>( geomV2 );
331  if ( !surface )
332  return 0.0;
333 
334  QgsPolygonV2* polygon = surface->surfaceToPolygon();
335 
336  double area = 0;
337  const QgsCurveV2* outerRing = polygon->exteriorRing();
338  area += measurePolygon( outerRing );
339 
340  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
341  {
342  const QgsCurveV2* innerRing = polygon->interiorRing( i );
343  area -= measurePolygon( innerRing );
344  }
345  delete polygon;
346  return area;
347  }
348  }
349 }
350 
351 double QgsDistanceArea::measure( const QgsGeometry *geometry ) const
352 {
353  if ( !geometry )
354  return 0.0;
355 
356  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
357  return measure( geomV2 );
358 }
359 
360 double QgsDistanceArea::measureArea( const QgsGeometry* geometry ) const
361 {
362  if ( !geometry )
363  return 0.0;
364 
365  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
366  return measure( geomV2, Area );
367 }
368 
369 double QgsDistanceArea::measureLength( const QgsGeometry* geometry ) const
370 {
371  if ( !geometry )
372  return 0.0;
373 
374  const QgsAbstractGeometryV2* 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  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
384  if ( !geomV2 || geomV2->dimension() < 2 )
385  {
386  return 0.0;
387  }
388 
389  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
390  {
391  return geomV2->perimeter();
392  }
393 
394  //create list with (single) surfaces
396  const QgsSurfaceV2* surf = dynamic_cast<const QgsSurfaceV2*>( geomV2 );
397  if ( surf )
398  {
399  surfaces.append( surf );
400  }
401  const QgsMultiSurfaceV2* multiSurf = dynamic_cast<const QgsMultiSurfaceV2*>( geomV2 );
402  if ( multiSurf )
403  {
404  surfaces.reserve(( surf ? 1 : 0 ) + multiSurf->numGeometries() );
405  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
406  {
407  surfaces.append( static_cast<const QgsSurfaceV2*>( multiSurf->geometryN( i ) ) );
408  }
409  }
410 
411  double length = 0;
413  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
414  {
415  if ( !*surfaceIt )
416  {
417  continue;
418  }
419 
420  QgsPolygonV2* poly = ( *surfaceIt )->surfaceToPolygon();
421  const QgsCurveV2* outerRing = poly->exteriorRing();
422  if ( outerRing )
423  {
424  length += measure( outerRing );
425  }
426  int nInnerRings = poly->numInteriorRings();
427  for ( int i = 0; i < nInnerRings; ++i )
428  {
429  length += measure( poly->interiorRing( i ) );
430  }
431  delete poly;
432  }
433  return length;
434 }
435 
436 double QgsDistanceArea::measureLine( const QgsCurveV2* curve ) const
437 {
438  if ( !curve )
439  {
440  return 0.0;
441  }
442 
443  QgsPointSequenceV2 linePointsV2;
444  QList<QgsPoint> linePoints;
445  curve->points( linePointsV2 );
446  QgsGeometry::convertPointList( linePointsV2, linePoints );
447  return measureLine( linePoints );
448 }
449 
450 double QgsDistanceArea::measureLine( const QList<QgsPoint> &points ) const
451 {
452  if ( points.size() < 2 )
453  return 0;
454 
455  double total = 0;
456  QgsPoint p1, p2;
457 
458  try
459  {
460  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
461  p1 = mCoordTransform->transform( points[0] );
462  else
463  p1 = points[0];
464 
465  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
466  {
467  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
468  {
469  p2 = mCoordTransform->transform( *i );
470  total += computeDistanceBearing( p1, p2 );
471  }
472  else
473  {
474  p2 = *i;
475  total += measureLine( p1, p2 );
476  }
477 
478  p1 = p2;
479  }
480 
481  return total;
482  }
483  catch ( QgsCsException &cse )
484  {
485  Q_UNUSED( cse );
486  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
487  return 0.0;
488  }
489 
490 }
491 
492 double QgsDistanceArea::measureLine( const QgsPoint &p1, const QgsPoint &p2 ) const
493 {
494  QGis::UnitType units;
495  return measureLine( p1, p2, units );
496 }
497 
498 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2, QGis::UnitType& units ) const
499 {
500  double result;
501  units = mCoordTransform->sourceCrs().mapUnits();
502 
503  try
504  {
505  QgsPoint pp1 = p1, pp2 = p2;
506 
507  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
508  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
509  {
510  units = QGis::Meters;
511  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
512  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform->sourceCrs().toProj4() ), 4 );
513  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform->destCRS().toProj4() ), 4 );
514  pp1 = mCoordTransform->transform( p1 );
515  pp2 = mCoordTransform->transform( p2 );
516  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
517  result = computeDistanceBearing( pp1, pp2 );
518  }
519  else
520  {
521  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
522  result = computeDistanceFlat( p1, p2 );
523  }
524  }
525  catch ( QgsCsException &cse )
526  {
527  Q_UNUSED( cse );
528  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
529  result = 0.0;
530  }
531  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
532  return result;
533 }
534 
536 {
537  return willUseEllipsoid() ? QGis::Meters : mCoordTransform->sourceCrs().mapUnits();
538 }
539 
541 {
543  QgsUnitTypes::distanceToAreaUnit( mCoordTransform->sourceCrs().mapUnits() );
544 }
545 
546 QgsConstWkbPtr QgsDistanceArea::measurePolygon( QgsConstWkbPtr wkbPtr, double* area, double* perimeter, bool hasZptr ) const
547 {
548  if ( !wkbPtr )
549  {
550  QgsDebugMsg( "no feature to measure" );
551  return wkbPtr;
552  }
553 
554  wkbPtr.readHeader();
555 
556  // get number of rings in the polygon
557  int numRings;
558  wkbPtr >> numRings;
559 
560  if ( numRings == 0 )
561  {
562  QgsDebugMsg( "no rings to measure" );
563  return QgsConstWkbPtr( nullptr, 0 );
564  }
565 
566  // Set pointer to the first ring
567  QList<QgsPoint> points;
568  QgsPoint pnt;
569  double x, y;
570  if ( area )
571  *area = 0;
572  if ( perimeter )
573  *perimeter = 0;
574 
575  try
576  {
577  for ( int idx = 0; idx < numRings; idx++ )
578  {
579  int nPoints;
580  wkbPtr >> nPoints;
581 
582  // Extract the points from the WKB and store in a pair of
583  // vectors.
584  for ( int jdx = 0; jdx < nPoints; jdx++ )
585  {
586  wkbPtr >> x >> y;
587  if ( hasZptr )
588  {
589  // totally ignore Z value
590  wkbPtr += sizeof( double );
591  }
592 
593  pnt = QgsPoint( x, y );
594 
595  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
596  {
597  pnt = mCoordTransform->transform( pnt );
598  }
599  points.append( pnt );
600  }
601 
602  if ( points.size() > 2 )
603  {
604  if ( area )
605  {
606  double areaTmp = computePolygonArea( points );
607  if ( idx == 0 )
608  {
609  // exterior ring
610  *area += areaTmp;
611  }
612  else
613  {
614  *area -= areaTmp; // interior rings
615  }
616  }
617 
618  if ( perimeter )
619  {
620  if ( idx == 0 )
621  {
622  // exterior ring
623  *perimeter += computeDistance( points );
624  }
625  }
626  }
627 
628  points.clear();
629  }
630  }
631  catch ( QgsCsException &cse )
632  {
633  Q_UNUSED( cse );
634  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
635  }
636 
637  return wkbPtr;
638 }
639 
640 double QgsDistanceArea::measurePolygon( const QgsCurveV2* curve ) const
641 {
642  if ( !curve )
643  {
644  return 0.0;
645  }
646 
647  QgsPointSequenceV2 linePointsV2;
648  curve->points( linePointsV2 );
649  QList<QgsPoint> linePoints;
650  QgsGeometry::convertPointList( linePointsV2, linePoints );
651  return measurePolygon( linePoints );
652 }
653 
654 
656 {
657  try
658  {
659  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
660  {
661  QList<QgsPoint> pts;
662  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
663  {
664  pts.append( mCoordTransform->transform( *i ) );
665  }
666  return computePolygonArea( pts );
667  }
668  else
669  {
670  return computePolygonArea( points );
671  }
672  }
673  catch ( QgsCsException &cse )
674  {
675  Q_UNUSED( cse );
676  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
677  return 0.0;
678  }
679 }
680 
681 
682 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 ) const
683 {
684  QgsPoint pp1 = p1, pp2 = p2;
685  double bearing;
686 
687  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
688  {
689  pp1 = mCoordTransform->transform( p1 );
690  pp2 = mCoordTransform->transform( p2 );
691  computeDistanceBearing( pp1, pp2, &bearing );
692  }
693  else //compute simple planar azimuth
694  {
695  double dx = p2.x() - p1.x();
696  double dy = p2.y() - p1.y();
697  bearing = atan2( dx, dy );
698  }
699 
700  return bearing;
701 }
702 
703 
705 // distance calculation
706 
708  const QgsPoint& p1, const QgsPoint& p2,
709  double* course1, double* course2 ) const
710 {
711  if ( qgsDoubleNear( p1.x(), p2.x() ) && qgsDoubleNear( p1.y(), p2.y() ) )
712  return 0;
713 
714  // ellipsoid
715  double a = mSemiMajor;
716  double b = mSemiMinor;
717  double f = 1 / mInvFlattening;
718 
719  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
720  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
721 
722  double L = p2_lon - p1_lon;
723  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
724  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
725  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
726  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
727  double lambda = L;
728  double lambdaP = 2 * M_PI;
729 
730  double sinLambda = 0;
731  double cosLambda = 0;
732  double sinSigma = 0;
733  double cosSigma = 0;
734  double sigma = 0;
735  double alpha = 0;
736  double cosSqAlpha = 0;
737  double cos2SigmaM = 0;
738  double C = 0;
739  double tu1 = 0;
740  double tu2 = 0;
741 
742  int iterLimit = 20;
743  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
744  {
745  sinLambda = sin( lambda );
746  cosLambda = cos( lambda );
747  tu1 = ( cosU2 * sinLambda );
748  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
749  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
750  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
751  sigma = atan2( sinSigma, cosSigma );
752  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
753  cosSqAlpha = cos( alpha ) * cos( alpha );
754  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
755  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
756  lambdaP = lambda;
757  lambda = L + ( 1 - C ) * f * sin( alpha ) *
758  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
759  }
760 
761  if ( iterLimit == 0 )
762  return -1; // formula failed to converge
763 
764  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
765  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
766  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
767  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
768  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
769  double s = b * A * ( sigma - deltaSigma );
770 
771  if ( course1 )
772  {
773  *course1 = atan2( tu1, tu2 );
774  }
775  if ( course2 )
776  {
777  // PI is added to return azimuth from P2 to P1
778  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
779  }
780 
781  return s;
782 }
783 
784 double QgsDistanceArea::computeDistanceFlat( const QgsPoint& p1, const QgsPoint& p2 ) const
785 {
786  return sqrt(( p2.x() - p1.x() ) * ( p2.x() - p1.x() ) + ( p2.y() - p1.y() ) * ( p2.y() - p1.y() ) );
787 }
788 
790 {
791  if ( points.size() < 2 )
792  return 0;
793 
794  double total = 0;
795  QgsPoint p1, p2;
796 
797  try
798  {
799  p1 = points[0];
800 
801  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
802  {
803  p2 = *i;
804  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
805  {
806  total += computeDistanceBearing( p1, p2 );
807  }
808  else
809  {
810  total += computeDistanceFlat( p1, p2 );
811  }
812 
813  p1 = p2;
814  }
815 
816  return total;
817  }
818  catch ( QgsCsException &cse )
819  {
820  Q_UNUSED( cse );
821  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
822  return 0.0;
823  }
824 }
825 
826 
827 
829 // stuff for measuring areas - copied from GRASS
830 // don't know how does it work, but it's working .)
831 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
832 
833 double QgsDistanceArea::getQ( double x ) const
834 {
835  double sinx, sinx2;
836 
837  sinx = sin( x );
838  sinx2 = sinx * sinx;
839 
840  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
841 }
842 
843 
844 double QgsDistanceArea::getQbar( double x ) const
845 {
846  double cosx, cosx2;
847 
848  cosx = cos( x );
849  cosx2 = cosx * cosx;
850 
851  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
852 }
853 
854 
856 {
857  //don't try to perform calculations if no ellipsoid
858  if ( mEllipsoid == GEO_NONE )
859  {
860  return;
861  }
862 
863  double a2 = ( mSemiMajor * mSemiMajor );
864  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
865  double e4, e6;
866 
867  m_TwoPI = M_PI + M_PI;
868 
869  e4 = e2 * e2;
870  e6 = e4 * e2;
871 
872  m_AE = a2 * ( 1 - e2 );
873 
874  m_QA = ( 2.0 / 3.0 ) * e2;
875  m_QB = ( 3.0 / 5.0 ) * e4;
876  m_QC = ( 4.0 / 7.0 ) * e6;
877 
878  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
879  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
880  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
881  m_QbarD = ( 4.0 / 49.0 ) * e6;
882 
883  m_Qp = getQ( M_PI / 2 );
884  m_E = 4 * M_PI * m_Qp * m_AE;
885  if ( m_E < 0.0 )
886  m_E = -m_E;
887 }
888 
889 
891 {
892  if ( points.isEmpty() )
893  {
894  return 0;
895  }
896 
897  double x1, y1, x2, y2, dx, dy;
898  double Qbar1, Qbar2;
899  double area;
900 
901  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
902  if (( ! mEllipsoidalMode ) || ( mEllipsoid == GEO_NONE ) )
903  {
904  return computePolygonFlatArea( points );
905  }
906  int n = points.size();
907  x2 = DEG2RAD( points[n-1].x() );
908  y2 = DEG2RAD( points[n-1].y() );
909  Qbar2 = getQbar( y2 );
910 
911  area = 0.0;
912 
913  for ( int i = 0; i < n; i++ )
914  {
915  x1 = x2;
916  y1 = y2;
917  Qbar1 = Qbar2;
918 
919  x2 = DEG2RAD( points[i].x() );
920  y2 = DEG2RAD( points[i].y() );
921  Qbar2 = getQbar( y2 );
922 
923  if ( x1 > x2 )
924  while ( x1 - x2 > M_PI )
925  x2 += m_TwoPI;
926  else if ( x2 > x1 )
927  while ( x2 - x1 > M_PI )
928  x1 += m_TwoPI;
929 
930  dx = x2 - x1;
931  area += dx * ( m_Qp - getQ( y2 ) );
932 
933  dy = y2 - y1;
934  if ( !qgsDoubleNear( dy, 0.0 ) )
935  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
936  }
937  if (( area *= m_AE ) < 0.0 )
938  area = -area;
939 
940  /* kludge - if polygon circles the south pole the area will be
941  * computed as if it cirlced the north pole. The correction is
942  * the difference between total surface area of the earth and
943  * the "north pole" area.
944  */
945  if ( area > m_E )
946  area = m_E;
947  if ( area > m_E / 2 )
948  area = m_E - area;
949 
950  return area;
951 }
952 
954 {
955  // Normal plane area calculations.
956  double area = 0.0;
957  int i, size;
958 
959  size = points.size();
960 
961  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
962  for ( i = 0; i < size; i++ )
963  {
964  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
965  // Using '% size', so that we always end with the starting point
966  // and thus close the polygon.
967  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
968  }
969  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
970  area = area / 2.0;
971  return qAbs( area ); // All areas are positive!
972 }
973 
974 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
975 {
976  QString unitLabel;
977 
978  switch ( u )
979  {
980  case QGis::Meters:
981  if ( isArea )
982  {
983  if ( keepBaseUnit )
984  {
985  unitLabel = QObject::trUtf8( " m²" );
986  }
987  else if ( qAbs( value ) > 1000000.0 )
988  {
989  unitLabel = QObject::trUtf8( " km²" );
990  value = value / 1000000.0;
991  }
992  else if ( qAbs( value ) > 10000.0 )
993  {
994  unitLabel = QObject::tr( " ha" );
995  value = value / 10000.0;
996  }
997  else
998  {
999  unitLabel = QObject::trUtf8( " m²" );
1000  }
1001  }
1002  else
1003  {
1004  if ( keepBaseUnit || qAbs( value ) == 0.0 )
1005  {
1006  unitLabel = QObject::tr( " m" );
1007  }
1008  else if ( qAbs( value ) > 1000.0 )
1009  {
1010  unitLabel = QObject::tr( " km" );
1011  value = value / 1000;
1012  }
1013  else if ( qAbs( value ) < 0.01 )
1014  {
1015  unitLabel = QObject::tr( " mm" );
1016  value = value * 1000;
1017  }
1018  else if ( qAbs( value ) < 0.1 )
1019  {
1020  unitLabel = QObject::tr( " cm" );
1021  value = value * 100;
1022  }
1023  else
1024  {
1025  unitLabel = QObject::tr( " m" );
1026  }
1027  }
1028  break;
1029  case QGis::Feet:
1030  if ( isArea )
1031  {
1032  if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
1033  {
1034  // < 0.5 acre show sq ft
1035  unitLabel = QObject::tr( " sq ft" );
1036  }
1037  else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
1038  {
1039  // < 0.5 sq mile show acre
1040  unitLabel = QObject::tr( " acres" );
1041  value /= 43560.0;
1042  }
1043  else
1044  {
1045  // above 0.5 acre show sq mi
1046  unitLabel = QObject::tr( " sq mile" );
1047  value /= 5280.0 * 5280.0;
1048  }
1049  }
1050  else
1051  {
1052  if ( qAbs( value ) <= 528.0 || keepBaseUnit )
1053  {
1054  if ( qAbs( value ) == 1.0 )
1055  {
1056  unitLabel = QObject::tr( " foot" );
1057  }
1058  else
1059  {
1060  unitLabel = QObject::tr( " feet" );
1061  }
1062  }
1063  else
1064  {
1065  unitLabel = QObject::tr( " mile" );
1066  value /= 5280.0;
1067  }
1068  }
1069  break;
1070  case QGis::NauticalMiles:
1071  if ( isArea )
1072  {
1073  unitLabel = QObject::tr( " sq. NM" );
1074  }
1075  else
1076  {
1077  unitLabel = QObject::tr( " NM" );
1078  }
1079  break;
1080  case QGis::Degrees:
1081  if ( isArea )
1082  {
1083  unitLabel = QObject::tr( " sq.deg." );
1084  }
1085  else
1086  {
1087  if ( qAbs( value ) == 1.0 )
1088  unitLabel = QObject::tr( " degree" );
1089  else
1090  unitLabel = QObject::tr( " degrees" );
1091  }
1092  break;
1093  case QGis::UnknownUnit:
1094  unitLabel.clear();
1095  break;
1096  default:
1097  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
1098  break;
1099  }
1100 
1101  return QString( "%L1%2" ).arg( value, 0, 'f', decimals ).arg( unitLabel );
1102 }
1103 
1104 QString QgsDistanceArea::formatDistance( double distance, int decimals, QGis::UnitType unit, bool keepBaseUnit )
1105 {
1106  QString unitLabel;
1107 
1108  switch ( unit )
1109  {
1110  case QGis::Meters:
1111  if ( keepBaseUnit || qAbs( distance ) == 0.0 )
1112  {
1113  unitLabel = QObject::tr( " m" );
1114  }
1115  else if ( qAbs( distance ) > 1000.0 )
1116  {
1117  unitLabel = QObject::tr( " km" );
1118  distance = distance / 1000;
1119  }
1120  else if ( qAbs( distance ) < 0.01 )
1121  {
1122  unitLabel = QObject::tr( " mm" );
1123  distance = distance * 1000;
1124  }
1125  else if ( qAbs( distance ) < 0.1 )
1126  {
1127  unitLabel = QObject::tr( " cm" );
1128  distance = distance * 100;
1129  }
1130  else
1131  {
1132  unitLabel = QObject::tr( " m" );
1133  }
1134  break;
1135 
1136  case QGis::Kilometers:
1137  if ( keepBaseUnit || qAbs( distance ) >= 1.0 )
1138  {
1139  unitLabel = QObject::tr( " km" );
1140  }
1141  else
1142  {
1143  unitLabel = QObject::tr( " m" );
1144  distance = distance * 1000;
1145  }
1146  break;
1147 
1148  case QGis::Feet:
1149  if ( qAbs( distance ) <= 5280.0 || keepBaseUnit )
1150  {
1151  unitLabel = QObject::tr( " ft" );
1152  }
1153  else
1154  {
1155  unitLabel = QObject::tr( " mi" );
1156  distance /= 5280.0;
1157  }
1158  break;
1159 
1160  case QGis::Yards:
1161  if ( qAbs( distance ) <= 1760.0 || keepBaseUnit )
1162  {
1163  unitLabel = QObject::tr( " yd" );
1164  }
1165  else
1166  {
1167  unitLabel = QObject::tr( " mi" );
1168  distance /= 1760.0;
1169  }
1170  break;
1171 
1172  case QGis::Miles:
1173  if ( qAbs( distance ) >= 1.0 || keepBaseUnit )
1174  {
1175  unitLabel = QObject::tr( " mi" );
1176  }
1177  else
1178  {
1179  unitLabel = QObject::tr( " ft" );
1180  distance *= 5280.0;
1181  }
1182  break;
1183 
1184  case QGis::NauticalMiles:
1185  unitLabel = QObject::tr( " NM" );
1186  break;
1187 
1188  case QGis::Degrees:
1189 
1190  if ( qAbs( distance ) == 1.0 )
1191  unitLabel = QObject::tr( " degree" );
1192  else
1193  unitLabel = QObject::tr( " degrees" );
1194  break;
1195 
1196  case QGis::UnknownUnit:
1197  unitLabel.clear();
1198  break;
1199  default:
1200  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( unit ) );
1201  break;
1202  }
1203 
1204  return QString( "%L1%2" ).arg( distance, 0, 'f', decimals ).arg( unitLabel );
1205 }
1206 
1207 QString QgsDistanceArea::formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit )
1208 {
1209  QString unitLabel;
1210 
1211  switch ( unit )
1212  {
1214  {
1215  if ( keepBaseUnit )
1216  {
1217  unitLabel = QObject::trUtf8( " m²" );
1218  }
1220  {
1221  unitLabel = QObject::trUtf8( " km²" );
1223  }
1225  {
1226  unitLabel = QObject::tr( " ha" );
1228  }
1229  else
1230  {
1231  unitLabel = QObject::trUtf8( " m²" );
1232  }
1233  break;
1234  }
1235 
1237  {
1238  unitLabel = QObject::trUtf8( " km²" );
1239  break;
1240  }
1241 
1243  {
1244  if ( keepBaseUnit )
1245  {
1246  unitLabel = QObject::trUtf8( " ft²" );
1247  }
1249  {
1250  unitLabel = QObject::trUtf8( " mi²" );
1252  }
1253  else
1254  {
1255  unitLabel = QObject::trUtf8( " ft²" );
1256  }
1257  break;
1258  }
1259 
1261  {
1262  if ( keepBaseUnit )
1263  {
1264  unitLabel = QObject::trUtf8( " yd²" );
1265  }
1267  {
1268  unitLabel = QObject::trUtf8( " mi²" );
1270  }
1271  else
1272  {
1273  unitLabel = QObject::trUtf8( " yd²" );
1274  }
1275  break;
1276  }
1277 
1279  {
1280  unitLabel = QObject::trUtf8( " mi²" );
1281  break;
1282  }
1283 
1285  {
1286  if ( keepBaseUnit )
1287  {
1288  unitLabel = QObject::trUtf8( " ha" );
1289  }
1291  {
1292  unitLabel = QObject::trUtf8( " km²" );
1294  }
1295  else
1296  {
1297  unitLabel = QObject::trUtf8( " ha" );
1298  }
1299  break;
1300  }
1301 
1302  case QgsUnitTypes::Acres:
1303  {
1304  if ( keepBaseUnit )
1305  {
1306  unitLabel = QObject::trUtf8( " ac" );
1307  }
1309  {
1310  unitLabel = QObject::trUtf8( " mi²" );
1312  }
1313  else
1314  {
1315  unitLabel = QObject::trUtf8( " ac" );
1316  }
1317  break;
1318  }
1319 
1321  {
1322  unitLabel = QObject::trUtf8( " nm²" );
1323  break;
1324  }
1325 
1327  {
1328  unitLabel = QObject::tr( " sq.deg." );
1329  break;
1330  }
1331 
1333  {
1334  unitLabel.clear();
1335  break;
1336  }
1337  }
1338 
1339  return QString( "%L1%2" ).arg( area, 0, 'f', decimals ).arg( unitLabel );
1340 }
1341 
1342 void QgsDistanceArea::convertMeasurement( double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea ) const
1343 {
1344  // Helper for converting between meters and feet and degrees and NauticalMiles...
1345  // The parameters measure and measureUnits are in/out
1346 
1347  if (( measureUnits == QGis::Degrees || measureUnits == QGis::Feet || measureUnits == QGis::NauticalMiles ) &&
1348  mEllipsoid != GEO_NONE &&
1349  mEllipsoidalMode )
1350  {
1351  // Measuring on an ellipsoid returned meters. Force!
1352  measureUnits = QGis::Meters;
1353  QgsDebugMsg( "We're measuring on an ellipsoid or using projections, the system is returning meters" );
1354  }
1355  else if ( mEllipsoidalMode && mEllipsoid == GEO_NONE )
1356  {
1357  // Measuring in plane within the source CRS. Force its map units
1358  measureUnits = mCoordTransform->sourceCrs().mapUnits();
1359  QgsDebugMsg( "We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1360  }
1361 
1362  // Gets the conversion factor between the specified units
1363  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, displayUnits );
1364  if ( isArea )
1365  factorUnits *= factorUnits;
1366 
1367  QgsDebugMsg( QString( "Converting %1 %2" ).arg( QString::number( measure ), QgsUnitTypes::toString( measureUnits ) ) );
1368  measure *= factorUnits;
1369  QgsDebugMsg( QString( "to %1 %2" ).arg( QString::number( measure ), QgsUnitTypes::toString( displayUnits ) ) );
1370  measureUnits = displayUnits;
1371 }
1372 
1373 double QgsDistanceArea::convertLengthMeasurement( double length, QGis::UnitType toUnits ) const
1374 {
1375  // get the conversion factor between the specified units
1376  QGis::UnitType measureUnits = lengthUnits();
1377  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1378 
1379  double result = length * factorUnits;
1380  QgsDebugMsg( QString( "Converted length of %1 %2 to %3 %4" ).arg( length )
1381  .arg( QgsUnitTypes::toString( measureUnits ) )
1382  .arg( result )
1383  .arg( QgsUnitTypes::toString( toUnits ) ) );
1384  return result;
1385 }
1386 
1388 {
1389  // get the conversion factor between the specified units
1390  QgsUnitTypes::AreaUnit measureUnits = areaUnits();
1391  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1392 
1393  double result = area * factorUnits;
1394  QgsDebugMsg( QString( "Converted area of %1 %2 to %3 %4" ).arg( area )
1395  .arg( QgsUnitTypes::toString( measureUnits ) )
1396  .arg( result )
1397  .arg( QgsUnitTypes::toString( toUnits ) ) );
1398  return result;
1399 }
const QgsCurveV2 * exteriorRing() const
QString ellipsoid() const
Returns ellipsoid&#39;s acronym.
QgsCoordinateReferenceSystem crsByProj4(const QString &proj4) const
Returns the CRS from a proj4 style formatted string.
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:352
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:38
Polygon geometry type.
Definition: qgspolygonv2.h:29
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)
QgsCoordinateReferenceSystem crsByOgcWmsCrs(const QString &ogcCrs) const
Returns the CRS from a given OGC WMS-format Coordinate Reference System string.
#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
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:460
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.
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 QgsLineStringV2 * 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...
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:159
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.
static QgsCRSCache * instance()
Returns a pointer to the QgsCRSCache singleton.
Definition: qgscrscache.cpp:91
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.
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.