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