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