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