QGIS API Documentation  2.99.0-Master (7d4f81d)
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 <QString>
18 #include <QObject>
19 
20 #include "qgsdistancearea.h"
21 #include "qgis.h"
22 #include "qgspointxy.h"
23 #include "qgscoordinatetransform.h"
25 #include "qgsgeometry.h"
26 #include "qgsgeometrycollection.h"
27 #include "qgslogger.h"
28 #include "qgsmessagelog.h"
29 #include "qgsmultisurface.h"
30 #include "qgswkbptr.h"
31 #include "qgslinestring.h"
32 #include "qgspolygon.h"
33 #include "qgssurface.h"
34 #include "qgsunittypes.h"
35 #include "qgsexception.h"
36 
37 // MSVC compiler doesn't have defined M_PI in math.h
38 #ifndef M_PI
39 #define M_PI 3.14159265358979323846
40 #endif
41 
42 #define DEG2RAD(x) ((x)*M_PI/180)
43 #define RAD2DEG(r) (180.0 * (r) / M_PI)
44 #define POW2(x) ((x)*(x))
45 
47 {
48  // init with default settings
49  mSemiMajor = -1.0;
50  mSemiMinor = -1.0;
51  mInvFlattening = -1.0;
54 }
55 
57 {
58  return mEllipsoid != GEO_NONE;
59 }
60 
62 {
63  mCoordTransform.setSourceCrs( srcCRS );
64 }
65 
67 {
68  // Shortcut if ellipsoid is none.
69  if ( ellipsoid == GEO_NONE )
70  {
71  mEllipsoid = GEO_NONE;
72  return true;
73  }
74 
76  if ( !params.valid )
77  {
78  return false;
79  }
80  else
81  {
82  mEllipsoid = ellipsoid;
83  setFromParams( params );
84  return true;
85  }
86 }
87 
88 // Inverse flattening is calculated with invf = a/(a-b)
89 // Also, b = a-(a/invf)
90 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
91 {
92  mEllipsoid = QStringLiteral( "PARAMETER:%1:%2" ).arg( qgsDoubleToString( semiMajor ) ).arg( qgsDoubleToString( semiMinor ) );
93  mSemiMajor = semiMajor;
94  mSemiMinor = semiMinor;
95  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
96 
97  computeAreaInit();
98 
99  return true;
100 }
101 
102 double QgsDistanceArea::measure( const QgsAbstractGeometry *geomV2, MeasureType type ) const
103 {
104  if ( !geomV2 )
105  {
106  return 0.0;
107  }
108 
109  int geomDimension = geomV2->dimension();
110  if ( geomDimension <= 0 )
111  {
112  return 0.0;
113  }
114 
115  MeasureType measureType = type;
116  if ( measureType == Default )
117  {
118  measureType = ( geomDimension == 1 ? Length : Area );
119  }
120 
121  if ( !willUseEllipsoid() )
122  {
123  //no transform required
124  if ( measureType == Length )
125  {
126  return geomV2->length();
127  }
128  else
129  {
130  return geomV2->area();
131  }
132  }
133  else
134  {
135  //multigeom is sum of measured parts
136  const QgsGeometryCollection *collection = qgsgeometry_cast<const QgsGeometryCollection *>( geomV2 );
137  if ( collection )
138  {
139  double sum = 0;
140  for ( int i = 0; i < collection->numGeometries(); ++i )
141  {
142  sum += measure( collection->geometryN( i ), measureType );
143  }
144  return sum;
145  }
146 
147  if ( measureType == Length )
148  {
149  const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( geomV2 );
150  if ( !curve )
151  {
152  return 0.0;
153  }
154 
155  QgsLineString *lineString = curve->curveToLine();
156  double length = measureLine( lineString );
157  delete lineString;
158  return length;
159  }
160  else
161  {
162  const QgsSurface *surface = qgsgeometry_cast<const QgsSurface *>( geomV2 );
163  if ( !surface )
164  return 0.0;
165 
166  QgsPolygonV2 *polygon = surface->surfaceToPolygon();
167 
168  double area = 0;
169  const QgsCurve *outerRing = polygon->exteriorRing();
170  area += measurePolygon( outerRing );
171 
172  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
173  {
174  const QgsCurve *innerRing = polygon->interiorRing( i );
175  area -= measurePolygon( innerRing );
176  }
177  delete polygon;
178  return area;
179  }
180  }
181 }
182 
183 double QgsDistanceArea::measureArea( const QgsGeometry &geometry ) const
184 {
185  if ( geometry.isNull() )
186  return 0.0;
187 
188  const QgsAbstractGeometry *geomV2 = geometry.geometry();
189  return measure( geomV2, Area );
190 }
191 
192 double QgsDistanceArea::measureLength( const QgsGeometry &geometry ) const
193 {
194  if ( geometry.isNull() )
195  return 0.0;
196 
197  const QgsAbstractGeometry *geomV2 = geometry.geometry();
198  return measure( geomV2, Length );
199 }
200 
201 double QgsDistanceArea::measurePerimeter( const QgsGeometry &geometry ) const
202 {
203  if ( geometry.isNull() )
204  return 0.0;
205 
206  const QgsAbstractGeometry *geomV2 = geometry.geometry();
207  if ( !geomV2 || geomV2->dimension() < 2 )
208  {
209  return 0.0;
210  }
211 
212  if ( !willUseEllipsoid() )
213  {
214  return geomV2->perimeter();
215  }
216 
217  //create list with (single) surfaces
218  QList< const QgsSurface * > surfaces;
219  const QgsSurface *surf = qgsgeometry_cast<const QgsSurface *>( geomV2 );
220  if ( surf )
221  {
222  surfaces.append( surf );
223  }
224  const QgsMultiSurface *multiSurf = qgsgeometry_cast<const QgsMultiSurface *>( geomV2 );
225  if ( multiSurf )
226  {
227  surfaces.reserve( ( surf ? 1 : 0 ) + multiSurf->numGeometries() );
228  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
229  {
230  surfaces.append( static_cast<const QgsSurface *>( multiSurf->geometryN( i ) ) );
231  }
232  }
233 
234  double length = 0;
235  QList<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
236  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
237  {
238  if ( !*surfaceIt )
239  {
240  continue;
241  }
242 
243  QgsPolygonV2 *poly = ( *surfaceIt )->surfaceToPolygon();
244  const QgsCurve *outerRing = poly->exteriorRing();
245  if ( outerRing )
246  {
247  length += measure( outerRing );
248  }
249  int nInnerRings = poly->numInteriorRings();
250  for ( int i = 0; i < nInnerRings; ++i )
251  {
252  length += measure( poly->interiorRing( i ) );
253  }
254  delete poly;
255  }
256  return length;
257 }
258 
259 double QgsDistanceArea::measureLine( const QgsCurve *curve ) const
260 {
261  if ( !curve )
262  {
263  return 0.0;
264  }
265 
266  QgsPointSequence linePointsV2;
267  QList<QgsPointXY> linePoints;
268  curve->points( linePointsV2 );
269  QgsGeometry::convertPointList( linePointsV2, linePoints );
270  return measureLine( linePoints );
271 }
272 
273 double QgsDistanceArea::measureLine( const QList<QgsPointXY> &points ) const
274 {
275  if ( points.size() < 2 )
276  return 0;
277 
278  double total = 0;
279  QgsPointXY p1, p2;
280 
281  try
282  {
283  if ( willUseEllipsoid() )
284  p1 = mCoordTransform.transform( points[0] );
285  else
286  p1 = points[0];
287 
288  for ( QList<QgsPointXY>::const_iterator i = points.begin(); i != points.end(); ++i )
289  {
290  if ( willUseEllipsoid() )
291  {
292  p2 = mCoordTransform.transform( *i );
293  total += computeDistanceBearing( p1, p2 );
294  }
295  else
296  {
297  p2 = *i;
298  total += measureLine( p1, p2 );
299  }
300 
301  p1 = p2;
302  }
303 
304  return total;
305  }
306  catch ( QgsCsException &cse )
307  {
308  Q_UNUSED( cse );
309  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
310  return 0.0;
311  }
312 
313 }
314 
315 double QgsDistanceArea::measureLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const
316 {
317  double result;
318 
319  try
320  {
321  QgsPointXY pp1 = p1, pp2 = p2;
322 
323  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
324  if ( willUseEllipsoid() )
325  {
326  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
327  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform.sourceCrs().toProj4() ), 4 );
328  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj4() ), 4 );
329  pp1 = mCoordTransform.transform( p1 );
330  pp2 = mCoordTransform.transform( p2 );
331  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
332  result = computeDistanceBearing( pp1, pp2 );
333  }
334  else
335  {
336  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
337  result = p2.distance( p1 );
338  }
339  }
340  catch ( QgsCsException &cse )
341  {
342  Q_UNUSED( cse );
343  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
344  result = 0.0;
345  }
346  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
347  return result;
348 }
349 
350 double QgsDistanceArea::measureLineProjected( const QgsPointXY &p1, double distance, double azimuth, QgsPointXY *projectedPoint ) const
351 {
352  double result = 0.0;
353  QgsPointXY p2;
354  if ( mCoordTransform.sourceCrs().isGeographic() && willUseEllipsoid() )
355  {
356  p2 = computeSpheroidProject( p1, distance, azimuth );
357  result = p1.distance( p2 );
358  }
359  else // Cartesian coordinates
360  {
361  result = distance; // Avoid rounding errors when using meters [return as sent]
362  if ( sourceCrs().mapUnits() != QgsUnitTypes::DistanceMeters )
363  {
364  distance = ( distance * QgsUnitTypes::fromUnitToUnitFactor( QgsUnitTypes::DistanceMeters, sourceCrs().mapUnits() ) );
365  result = p1.distance( p2 );
366  }
367  p2 = p1.project( distance, azimuth );
368  }
369  QgsDebugMsgLevel( QString( "Converted distance of %1 %2 to %3 distance %4 %5, using azimuth[%6] from point[%7] to point[%8] sourceCrs[%9] mEllipsoid[%10] isGeographic[%11] [%12]" )
370  .arg( QString::number( distance, 'f', 7 ) )
372  .arg( QString::number( result, 'f', 7 ) )
373  .arg( ( ( mCoordTransform.sourceCrs().isGeographic() ) == 1 ? QString( "Geographic" ) : QString( "Cartesian" ) ) )
374  .arg( QgsUnitTypes::toString( sourceCrs().mapUnits() ) )
375  .arg( azimuth )
376  .arg( p1.wellKnownText() )
377  .arg( p2.wellKnownText() )
378  .arg( sourceCrs().description() )
379  .arg( mEllipsoid )
380  .arg( sourceCrs().isGeographic() )
381  .arg( QString( "SemiMajor[%1] SemiMinor[%2] InvFlattening[%3] " ).arg( QString::number( mSemiMajor, 'f', 7 ) ).arg( QString::number( mSemiMinor, 'f', 7 ) ).arg( QString::number( mInvFlattening, 'f', 7 ) ) ), 4 );
382  if ( projectedPoint )
383  {
384  *projectedPoint = QgsPointXY( p2 );
385  }
386  return result;
387 }
388 
389 /*
390  * From original rttopo documentation:
391  * Tested against:
392  * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
393  * and
394  * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
395  */
397  const QgsPointXY &p1, double distance, double azimuth ) const
398 {
399  // ellipsoid
400  double a = mSemiMajor;
401  double b = mSemiMinor;
402  double f = 1 / mInvFlattening;
403  if ( ( ( a < 0 ) && ( b < 0 ) ) ||
404  ( ( p1.x() < -180.0 ) || ( p1.x() > 180.0 ) || ( p1.y() < -85.05115 ) || ( p1.y() > 85.05115 ) ) )
405  {
406  // latitudes outside these bounds cause the calculations to become unstable and can return invalid results
407  return QgsPoint( 0, 0 );
408 
409  }
410  double radians_lat = DEG2RAD( p1.y() );
411  double radians_long = DEG2RAD( p1.x() );
412  double b2 = POW2( b ); // spheroid_mu2
413  double omf = 1 - f;
414  double tan_u1 = omf * tan( radians_lat );
415  double u1 = atan( tan_u1 );
416  double sigma, last_sigma, delta_sigma, two_sigma_m;
417  double sigma1, sin_alpha, alpha, cos_alphasq;
418  double u2, A, B;
419  double lat2, lambda, lambda2, C, omega;
420  int i = 0;
421  if ( azimuth < 0.0 )
422  {
423  azimuth = azimuth + M_PI * 2.0;
424  }
425  if ( azimuth > ( M_PI * 2.0 ) )
426  {
427  azimuth = azimuth - M_PI * 2.0;
428  }
429  sigma1 = atan2( tan_u1, cos( azimuth ) );
430  sin_alpha = cos( u1 ) * sin( azimuth );
431  alpha = asin( sin_alpha );
432  cos_alphasq = 1.0 - POW2( sin_alpha );
433  u2 = POW2( cos( alpha ) ) * ( POW2( a ) - b2 ) / b2; // spheroid_mu2
434  A = 1.0 + ( u2 / 16384.0 ) * ( 4096.0 + u2 * ( -768.0 + u2 * ( 320.0 - 175.0 * u2 ) ) );
435  B = ( u2 / 1024.0 ) * ( 256.0 + u2 * ( -128.0 + u2 * ( 74.0 - 47.0 * u2 ) ) );
436  sigma = ( distance / ( b * A ) );
437  do
438  {
439  two_sigma_m = 2.0 * sigma1 + sigma;
440  delta_sigma = B * sin( sigma ) * ( cos( two_sigma_m ) + ( B / 4.0 ) * ( cos( sigma ) * ( -1.0 + 2.0 * POW2( cos( two_sigma_m ) ) - ( B / 6.0 ) * cos( two_sigma_m ) * ( -3.0 + 4.0 * POW2( sin( sigma ) ) ) * ( -3.0 + 4.0 * POW2( cos( two_sigma_m ) ) ) ) ) );
441  last_sigma = sigma;
442  sigma = ( distance / ( b * A ) ) + delta_sigma;
443  i++;
444  }
445  while ( i < 999 && qAbs( ( last_sigma - sigma ) / sigma ) > 1.0e-9 );
446 
447  lat2 = atan2( ( sin( u1 ) * cos( sigma ) + cos( u1 ) * sin( sigma ) *
448  cos( azimuth ) ), ( omf * sqrt( POW2( sin_alpha ) +
449  POW2( sin( u1 ) * sin( sigma ) - cos( u1 ) * cos( sigma ) *
450  cos( azimuth ) ) ) ) );
451  lambda = atan2( ( sin( sigma ) * sin( azimuth ) ), ( cos( u1 ) * cos( sigma ) -
452  sin( u1 ) * sin( sigma ) * cos( azimuth ) ) );
453  C = ( f / 16.0 ) * cos_alphasq * ( 4.0 + f * ( 4.0 - 3.0 * cos_alphasq ) );
454  omega = lambda - ( 1.0 - C ) * f * sin_alpha * ( sigma + C * sin( sigma ) *
455  ( cos( two_sigma_m ) + C * cos( sigma ) * ( -1.0 + 2.0 * POW2( cos( two_sigma_m ) ) ) ) );
456  lambda2 = radians_long + omega;
457  return QgsPointXY( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
458 }
459 
461 {
462  return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();
463 }
464 
466 {
468  QgsUnitTypes::distanceToAreaUnit( mCoordTransform.sourceCrs().mapUnits() );
469 }
470 
471 double QgsDistanceArea::measurePolygon( const QgsCurve *curve ) const
472 {
473  if ( !curve )
474  {
475  return 0.0;
476  }
477 
478  QgsPointSequence linePointsV2;
479  curve->points( linePointsV2 );
480  QList<QgsPointXY> linePoints;
481  QgsGeometry::convertPointList( linePointsV2, linePoints );
482  return measurePolygon( linePoints );
483 }
484 
485 
486 double QgsDistanceArea::measurePolygon( const QList<QgsPointXY> &points ) const
487 {
488  try
489  {
490  if ( willUseEllipsoid() )
491  {
492  QList<QgsPointXY> pts;
493  for ( QList<QgsPointXY>::const_iterator i = points.begin(); i != points.end(); ++i )
494  {
495  pts.append( mCoordTransform.transform( *i ) );
496  }
497  return computePolygonArea( pts );
498  }
499  else
500  {
501  return computePolygonArea( points );
502  }
503  }
504  catch ( QgsCsException &cse )
505  {
506  Q_UNUSED( cse );
507  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
508  return 0.0;
509  }
510 }
511 
512 
513 double QgsDistanceArea::bearing( const QgsPointXY &p1, const QgsPointXY &p2 ) const
514 {
515  QgsPointXY pp1 = p1, pp2 = p2;
516  double bearing;
517 
518  if ( willUseEllipsoid() )
519  {
520  pp1 = mCoordTransform.transform( p1 );
521  pp2 = mCoordTransform.transform( p2 );
522  computeDistanceBearing( pp1, pp2, &bearing );
523  }
524  else //compute simple planar azimuth
525  {
526  double dx = p2.x() - p1.x();
527  double dy = p2.y() - p1.y();
528  bearing = atan2( dx, dy );
529  }
530 
531  return bearing;
532 }
533 
534 
536 // distance calculation
537 
538 double QgsDistanceArea::computeDistanceBearing(
539  const QgsPointXY &p1, const QgsPointXY &p2,
540  double *course1, double *course2 ) const
541 {
542  if ( qgsDoubleNear( p1.x(), p2.x() ) && qgsDoubleNear( p1.y(), p2.y() ) )
543  return 0;
544 
545  // ellipsoid
546  double a = mSemiMajor;
547  double b = mSemiMinor;
548  double f = 1 / mInvFlattening;
549 
550  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
551  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
552 
553  double L = p2_lon - p1_lon;
554  double U1 = atan( ( 1 - f ) * tan( p1_lat ) );
555  double U2 = atan( ( 1 - f ) * tan( p2_lat ) );
556  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
557  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
558  double lambda = L;
559  double lambdaP = 2 * M_PI;
560 
561  double sinLambda = 0;
562  double cosLambda = 0;
563  double sinSigma = 0;
564  double cosSigma = 0;
565  double sigma = 0;
566  double alpha = 0;
567  double cosSqAlpha = 0;
568  double cos2SigmaM = 0;
569  double C = 0;
570  double tu1 = 0;
571  double tu2 = 0;
572 
573  int iterLimit = 20;
574  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
575  {
576  sinLambda = sin( lambda );
577  cosLambda = cos( lambda );
578  tu1 = ( cosU2 * sinLambda );
579  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
580  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
581  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
582  sigma = atan2( sinSigma, cosSigma );
583  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
584  cosSqAlpha = cos( alpha ) * cos( alpha );
585  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
586  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
587  lambdaP = lambda;
588  lambda = L + ( 1 - C ) * f * sin( alpha ) *
589  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
590  }
591 
592  if ( iterLimit == 0 )
593  return -1; // formula failed to converge
594 
595  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
596  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
597  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
598  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
599  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
600  double s = b * A * ( sigma - deltaSigma );
601 
602  if ( course1 )
603  {
604  *course1 = atan2( tu1, tu2 );
605  }
606  if ( course2 )
607  {
608  // PI is added to return azimuth from P2 to P1
609  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
610  }
611 
612  return s;
613 }
614 
616 // stuff for measuring areas - copied from GRASS
617 // don't know how does it work, but it's working .)
618 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
619 
620 double QgsDistanceArea::getQ( double x ) const
621 {
622  double sinx, sinx2;
623 
624  sinx = sin( x );
625  sinx2 = sinx * sinx;
626 
627  return sinx * ( 1 + sinx2 * ( m_QA + sinx2 * ( m_QB + sinx2 * m_QC ) ) );
628 }
629 
630 
631 double QgsDistanceArea::getQbar( double x ) const
632 {
633  double cosx, cosx2;
634 
635  cosx = cos( x );
636  cosx2 = cosx * cosx;
637 
638  return cosx * ( m_QbarA + cosx2 * ( m_QbarB + cosx2 * ( m_QbarC + cosx2 * m_QbarD ) ) );
639 }
640 
641 
642 void QgsDistanceArea::computeAreaInit()
643 {
644  //don't try to perform calculations if no ellipsoid
645  if ( mEllipsoid == GEO_NONE )
646  {
647  return;
648  }
649 
650  double a2 = ( mSemiMajor * mSemiMajor );
651  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
652  double e4, e6;
653 
654  m_TwoPI = M_PI + M_PI;
655 
656  e4 = e2 * e2;
657  e6 = e4 * e2;
658 
659  m_AE = a2 * ( 1 - e2 );
660 
661  m_QA = ( 2.0 / 3.0 ) * e2;
662  m_QB = ( 3.0 / 5.0 ) * e4;
663  m_QC = ( 4.0 / 7.0 ) * e6;
664 
665  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
666  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
667  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
668  m_QbarD = ( 4.0 / 49.0 ) * e6;
669 
670  m_Qp = getQ( M_PI / 2 );
671  m_E = 4 * M_PI * m_Qp * m_AE;
672  if ( m_E < 0.0 )
673  m_E = -m_E;
674 }
675 
676 void QgsDistanceArea::setFromParams( const QgsEllipsoidUtils::EllipsoidParameters &params )
677 {
678  if ( params.useCustomParameters )
679  {
680  setEllipsoid( params.semiMajor, params.semiMinor );
681  }
682  else
683  {
684  mSemiMajor = params.semiMajor;
685  mSemiMinor = params.semiMinor;
686  mInvFlattening = params.inverseFlattening;
687  mCoordTransform.setDestinationCrs( params.crs );
688  // precalculate some values for area calculations
689  computeAreaInit();
690  }
691 }
692 
693 double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) const
694 {
695  if ( points.isEmpty() )
696  {
697  return 0;
698  }
699 
700  // IMPORTANT
701  // don't change anything here without reporting the changes to upstream (GRASS)
702  // let's all be good opensource citizens and share the improvements!
703 
704  double x1, y1, x2, y2, dx, dy;
705  double Qbar1, Qbar2;
706  double area;
707 
708  /* GRASS comment: threshold for dy, should be between 1e-4 and 1e-7
709  * See relevant discussion at https://trac.osgeo.org/grass/ticket/3369
710  */
711  const double thresh = 1e-6;
712 
713  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
714  if ( !willUseEllipsoid() )
715  {
716  return computePolygonFlatArea( points );
717  }
718  int n = points.size();
719  x2 = DEG2RAD( points[n - 1].x() );
720  y2 = DEG2RAD( points[n - 1].y() );
721  Qbar2 = getQbar( y2 );
722 
723  area = 0.0;
724 
725  for ( int i = 0; i < n; i++ )
726  {
727  x1 = x2;
728  y1 = y2;
729  Qbar1 = Qbar2;
730 
731  x2 = DEG2RAD( points[i].x() );
732  y2 = DEG2RAD( points[i].y() );
733  Qbar2 = getQbar( y2 );
734 
735  if ( x1 > x2 )
736  while ( x1 - x2 > M_PI )
737  x2 += m_TwoPI;
738  else if ( x2 > x1 )
739  while ( x2 - x1 > M_PI )
740  x1 += m_TwoPI;
741 
742  dx = x2 - x1;
743  dy = y2 - y1;
744  if ( qAbs( dy ) > thresh )
745  {
746  /* account for different latitudes y1, y2 */
747  area += dx * ( m_Qp - ( Qbar2 - Qbar1 ) / dy );
748  }
749  else
750  {
751  /* latitudes y1, y2 are (nearly) identical */
752 
753  /* if y2 becomes similar to y1, i.e. y2 -> y1
754  * Qbar2 - Qbar1 -> 0 and dy -> 0
755  * (Qbar2 - Qbar1) / dy -> ?
756  * (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
757  * Metz 2017
758  */
759  area += dx * ( m_Qp - getQ( ( y1 + y2 ) / 2.0 ) );
760 
761  /* original:
762  * area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
763  */
764  }
765  }
766  if ( ( area *= m_AE ) < 0.0 )
767  area = -area;
768 
769  /* kludge - if polygon circles the south pole the area will be
770  * computed as if it cirlced the north pole. The correction is
771  * the difference between total surface area of the earth and
772  * the "north pole" area.
773  */
774  if ( area > m_E )
775  area = m_E;
776  if ( area > m_E / 2 )
777  area = m_E - area;
778 
779  return area;
780 }
781 
782 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPointXY> &points ) const
783 {
784  // Normal plane area calculations.
785  double area = 0.0;
786  int i, size;
787 
788  size = points.size();
789 
790  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
791  for ( i = 0; i < size; i++ )
792  {
793  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
794  // Using '% size', so that we always end with the starting point
795  // and thus close the polygon.
796  area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
797  }
798  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
799  area = area / 2.0;
800  return qAbs( area ); // All areas are positive!
801 }
802 
803 QString QgsDistanceArea::formatDistance( double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit )
804 {
805  return QgsUnitTypes::formatDistance( distance, decimals, unit, keepBaseUnit );
806 }
807 
808 QString QgsDistanceArea::formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit )
809 {
810  return QgsUnitTypes::formatArea( area, decimals, unit, keepBaseUnit );
811 }
812 
814 {
815  // get the conversion factor between the specified units
816  QgsUnitTypes::DistanceUnit measureUnits = lengthUnits();
817  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
818 
819  double result = length * factorUnits;
820  QgsDebugMsgLevel( QString( "Converted length of %1 %2 to %3 %4" ).arg( length )
821  .arg( QgsUnitTypes::toString( measureUnits ) )
822  .arg( result )
823  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
824  return result;
825 }
826 
828 {
829  // get the conversion factor between the specified units
830  QgsUnitTypes::AreaUnit measureUnits = areaUnits();
831  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
832 
833  double result = area * factorUnits;
834  QgsDebugMsgLevel( QString( "Converted area of %1 %2 to %3 %4" ).arg( area )
835  .arg( QgsUnitTypes::toString( measureUnits ) )
836  .arg( result )
837  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
838  return result;
839 }
const QgsCurve * interiorRing(int i) const
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)
static Q_INVOKABLE QString toString(QgsUnitTypes::DistanceUnit unit)
Returns a translated string representing a distance unit.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
virtual QgsPolygonV2 * surfaceToPolygon() const =0
Get a polygon representation of this surface.
QString wellKnownText() const
Return the well known text representation for the point.
Definition: qgspointxy.cpp:259
double y
Definition: qgspointxy.h:47
A class to represent a 2D point.
Definition: qgspointxy.h:42
double measurePolygon(const QList< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
#define DEG2RAD(x)
QString toProj4() const
Returns a Proj4 string representation of this CRS.
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:96
double convertLengthMeasurement(double length, QgsUnitTypes::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
double convertAreaMeasurement(double area, QgsUnitTypes::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
Contains parameters for an ellipsoid.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
Definition: qgis.h:210
Multi surface geometry collection.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:71
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
bool willUseEllipsoid() const
Returns whether calculations will use the ellipsoid.
Polygon geometry type.
Definition: qgspolygon.h:30
static Q_INVOKABLE QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
int numInteriorRings() const
virtual double length() const
Returns the length of the geometry.
QgsPointXY computeSpheroidProject(const QgsPointXY &p1, double distance=1, double azimuth=M_PI/2) const
Given a location, an azimuth and a distance, computes the location of the projected point...
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from...
static void convertPointList(const QList< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPointV2.
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
bool valid
Whether ellipsoid parameters are valid.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:38
Geometry collection.
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:198
#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:323
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
QString description() const
Returns the descriptive name of the CRS, e.g., "WGS 84" or "GDA 94 / Vicgrid94".
QgsPointXY project(double distance, double bearing) const
Returns a new point which corresponds to this point projected by a specified distance in a specified ...
Definition: qgspointxy.cpp:291
virtual double area() const
Returns the area of the geometry.
bool isGeographic() const
Returns whether the CRS is a geographic CRS (using lat/lon coordinates)
Abstract base class for curved geometry type.
Definition: qgscurve.h:34
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
Abstract base class for all geometries.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
const QgsCurve * exteriorRing() const
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.cpp:274
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:36
#define POW2(x)
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
double x
Definition: qgspointxy.h:46
double measureLineProjected(const QgsPointXY &p1, double distance=1, double azimuth=M_PI/2, QgsPointXY *projectedPoint=nullptr) const
Calculates the distance from one point with distance in meters and azimuth (direction) When the sourc...
int numGeometries() const
Returns the number of geometries within the collection.
DistanceUnit
Units of distance.
Definition: qgsunittypes.h:42
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
QgsUnitTypes::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
static void logMessage(const QString &message, const QString &tag=QString(), MessageLevel level=QgsMessageLog::WARNING)
add a message to the instance (and create it if necessary)
QString ellipsoid() const
Returns ellipsoid&#39;s acronym.
QString toString() const
String representation of the point (x,y)
Definition: qgspointxy.cpp:45
static Q_INVOKABLE QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const
Transform the point from the source CRS to the destination CRS.
QgsPolygonV2 * surfaceToPolygon() const override
Get a polygon representation of this surface.
Definition: qgspolygon.cpp:298
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:40
virtual double perimeter() const
Returns the perimeter of the geometry.
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
This class represents a coordinate reference system (CRS).
double inverseFlattening
Inverse flattening.
static QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:64
void setSourceCrs(const QgsCoordinateReferenceSystem &srcCRS)
Sets source spatial reference system.
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
QList< QgsPoint > QgsPointSequence
static Q_INVOKABLE double fromUnitToUnitFactor(QgsUnitTypes::DistanceUnit fromUnit, QgsUnitTypes::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
QgsDistanceArea()
Constructor.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
double measureLine(const QList< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
static Q_INVOKABLE AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
QgsAbstractGeometry * geometry() const
Returns the underlying geometry store.
AreaUnit
Units of area.
Definition: qgsunittypes.h:66
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
static QgsCoordinateReferenceSystem fromSrsId(long srsId)
Creates a CRS from a specified QGIS SRS ID.
#define RAD2DEG(r)
static QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.