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