QGIS API Documentation  2.99.0-Master (d55fa22)
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 "qgspoint.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 "qgscsexception.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
51 }
52 
54 {
55  return mEllipsoid != GEO_NONE;
56 }
57 
59 {
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 ) ).arg( 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 = dynamic_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 = dynamic_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 = dynamic_cast<const QgsSurface *>( geomV2 );
160  if ( !surface )
161  return 0.0;
162 
163  QgsPolygonV2 *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.geometry();
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.geometry();
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.geometry();
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  QList< const QgsSurface * > surfaces;
216  const QgsSurface *surf = dynamic_cast<const QgsSurface *>( geomV2 );
217  if ( surf )
218  {
219  surfaces.append( surf );
220  }
221  const QgsMultiSurface *multiSurf = dynamic_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  QList<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
233  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
234  {
235  if ( !*surfaceIt )
236  {
237  continue;
238  }
239 
240  QgsPolygonV2 *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  QList<QgsPoint> linePoints;
265  curve->points( linePointsV2 );
266  QgsGeometry::convertPointList( linePointsV2, linePoints );
267  return measureLine( linePoints );
268 }
269 
270 double QgsDistanceArea::measureLine( const QList<QgsPoint> &points ) const
271 {
272  if ( points.size() < 2 )
273  return 0;
274 
275  double total = 0;
276  QgsPoint p1, p2;
277 
278  try
279  {
280  if ( willUseEllipsoid() )
281  p1 = mCoordTransform.transform( points[0] );
282  else
283  p1 = points[0];
284 
285  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++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 QgsPoint &p1, const QgsPoint &p2 ) const
313 {
314  double result;
315 
316  try
317  {
318  QgsPoint pp1 = p1, pp2 = p2;
319 
320  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
321  if ( willUseEllipsoid() )
322  {
323  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
324  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform.sourceCrs().toProj4() ), 4 );
325  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj4() ), 4 );
326  pp1 = mCoordTransform.transform( p1 );
327  pp2 = mCoordTransform.transform( p2 );
328  QgsDebugMsgLevel( QString( "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( "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( QString( "The result was %1" ).arg( result ), 3 );
344  return result;
345 }
346 
347 double QgsDistanceArea::measureLineProjected( const QgsPoint &p1, double distance, double azimuth, QgsPoint *projectedPoint ) const
348 {
349  double result = 0.0;
350  QgsPoint 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  if ( projectedPoint )
367  {
368  *projectedPoint = QgsPoint( p2 );
369  }
370  return result;
371 }
372 
373 /*
374  * From original rttopo documentation:
375  * Tested against:
376  * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
377  * and
378  * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
379  */
381  const QgsPoint &p1, double distance, double azimuth ) const
382 {
383  // ellipsoid
384  double a = mSemiMajor;
385  double b = mSemiMinor;
386  double f = 1 / mInvFlattening;
387  double radians_lat = DEG2RAD( p1.y() );
388  double radians_long = DEG2RAD( p1.x() );
389  double b2 = POW2( b ); // spheroid_mu2
390  double omf = 1 - f;
391  double tan_u1 = omf * tan( radians_lat );
392  double u1 = atan( tan_u1 );
393  double sigma, last_sigma, delta_sigma, two_sigma_m;
394  double sigma1, sin_alpha, alpha, cos_alphasq;
395  double u2, A, B;
396  double lat2, lambda, lambda2, C, omega;
397  int i = 0;
398  if ( azimuth < 0.0 )
399  {
400  azimuth = azimuth + M_PI * 2.0;
401  }
402  if ( azimuth > ( M_PI * 2.0 ) )
403  {
404  azimuth = azimuth - M_PI * 2.0;
405  }
406  sigma1 = atan2( tan_u1, cos( azimuth ) );
407  sin_alpha = cos( u1 ) * sin( azimuth );
408  alpha = asin( sin_alpha );
409  cos_alphasq = 1.0 - POW2( sin_alpha );
410  u2 = POW2( cos( alpha ) ) * ( POW2( a ) - b2 ) / b2; // spheroid_mu2
411  A = 1.0 + ( u2 / 16384.0 ) * ( 4096.0 + u2 * ( -768.0 + u2 * ( 320.0 - 175.0 * u2 ) ) );
412  B = ( u2 / 1024.0 ) * ( 256.0 + u2 * ( -128.0 + u2 * ( 74.0 - 47.0 * u2 ) ) );
413  sigma = ( distance / ( b * A ) );
414  do
415  {
416  two_sigma_m = 2.0 * sigma1 + sigma;
417  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 ) ) ) ) ) );
418  last_sigma = sigma;
419  sigma = ( distance / ( b * A ) ) + delta_sigma;
420  i++;
421  }
422  while ( i < 999 && qAbs( ( last_sigma - sigma ) / sigma ) > 1.0e-9 );
423 
424  lat2 = atan2( ( sin( u1 ) * cos( sigma ) + cos( u1 ) * sin( sigma ) *
425  cos( azimuth ) ), ( omf * sqrt( POW2( sin_alpha ) +
426  POW2( sin( u1 ) * sin( sigma ) - cos( u1 ) * cos( sigma ) *
427  cos( azimuth ) ) ) ) );
428  lambda = atan2( ( sin( sigma ) * sin( azimuth ) ), ( cos( u1 ) * cos( sigma ) -
429  sin( u1 ) * sin( sigma ) * cos( azimuth ) ) );
430  C = ( f / 16.0 ) * cos_alphasq * ( 4.0 + f * ( 4.0 - 3.0 * cos_alphasq ) );
431  omega = lambda - ( 1.0 - C ) * f * sin_alpha * ( sigma + C * sin( sigma ) *
432  ( cos( two_sigma_m ) + C * cos( sigma ) * ( -1.0 + 2.0 * POW2( cos( two_sigma_m ) ) ) ) );
433  lambda2 = radians_long + omega;
434  return QgsPoint( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
435 }
436 
438 {
439  return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();
440 }
441 
443 {
445  QgsUnitTypes::distanceToAreaUnit( mCoordTransform.sourceCrs().mapUnits() );
446 }
447 
448 double QgsDistanceArea::measurePolygon( const QgsCurve *curve ) const
449 {
450  if ( !curve )
451  {
452  return 0.0;
453  }
454 
455  QgsPointSequence linePointsV2;
456  curve->points( linePointsV2 );
457  QList<QgsPoint> linePoints;
458  QgsGeometry::convertPointList( linePointsV2, linePoints );
459  return measurePolygon( linePoints );
460 }
461 
462 
463 double QgsDistanceArea::measurePolygon( const QList<QgsPoint> &points ) const
464 {
465  try
466  {
467  if ( willUseEllipsoid() )
468  {
469  QList<QgsPoint> pts;
470  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
471  {
472  pts.append( mCoordTransform.transform( *i ) );
473  }
474  return computePolygonArea( pts );
475  }
476  else
477  {
478  return computePolygonArea( points );
479  }
480  }
481  catch ( QgsCsException &cse )
482  {
483  Q_UNUSED( cse );
484  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
485  return 0.0;
486  }
487 }
488 
489 
490 double QgsDistanceArea::bearing( const QgsPoint &p1, const QgsPoint &p2 ) const
491 {
492  QgsPoint pp1 = p1, pp2 = p2;
493  double bearing;
494 
495  if ( willUseEllipsoid() )
496  {
497  pp1 = mCoordTransform.transform( p1 );
498  pp2 = mCoordTransform.transform( p2 );
499  computeDistanceBearing( pp1, pp2, &bearing );
500  }
501  else //compute simple planar azimuth
502  {
503  double dx = p2.x() - p1.x();
504  double dy = p2.y() - p1.y();
505  bearing = atan2( dx, dy );
506  }
507 
508  return bearing;
509 }
510 
511 
513 // distance calculation
514 
515 double QgsDistanceArea::computeDistanceBearing(
516  const QgsPoint &p1, const QgsPoint &p2,
517  double *course1, double *course2 ) const
518 {
519  if ( qgsDoubleNear( p1.x(), p2.x() ) && qgsDoubleNear( p1.y(), p2.y() ) )
520  return 0;
521 
522  // ellipsoid
523  double a = mSemiMajor;
524  double b = mSemiMinor;
525  double f = 1 / mInvFlattening;
526 
527  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
528  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
529 
530  double L = p2_lon - p1_lon;
531  double U1 = atan( ( 1 - f ) * tan( p1_lat ) );
532  double U2 = atan( ( 1 - f ) * tan( p2_lat ) );
533  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
534  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
535  double lambda = L;
536  double lambdaP = 2 * M_PI;
537 
538  double sinLambda = 0;
539  double cosLambda = 0;
540  double sinSigma = 0;
541  double cosSigma = 0;
542  double sigma = 0;
543  double alpha = 0;
544  double cosSqAlpha = 0;
545  double cos2SigmaM = 0;
546  double C = 0;
547  double tu1 = 0;
548  double tu2 = 0;
549 
550  int iterLimit = 20;
551  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
552  {
553  sinLambda = sin( lambda );
554  cosLambda = cos( lambda );
555  tu1 = ( cosU2 * sinLambda );
556  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
557  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
558  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
559  sigma = atan2( sinSigma, cosSigma );
560  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
561  cosSqAlpha = cos( alpha ) * cos( alpha );
562  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
563  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
564  lambdaP = lambda;
565  lambda = L + ( 1 - C ) * f * sin( alpha ) *
566  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
567  }
568 
569  if ( iterLimit == 0 )
570  return -1; // formula failed to converge
571 
572  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
573  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
574  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
575  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
576  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
577  double s = b * A * ( sigma - deltaSigma );
578 
579  if ( course1 )
580  {
581  *course1 = atan2( tu1, tu2 );
582  }
583  if ( course2 )
584  {
585  // PI is added to return azimuth from P2 to P1
586  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
587  }
588 
589  return s;
590 }
591 
593 // stuff for measuring areas - copied from GRASS
594 // don't know how does it work, but it's working .)
595 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
596 
597 double QgsDistanceArea::getQ( double x ) const
598 {
599  double sinx, sinx2;
600 
601  sinx = sin( x );
602  sinx2 = sinx * sinx;
603 
604  return sinx * ( 1 + sinx2 * ( m_QA + sinx2 * ( m_QB + sinx2 * m_QC ) ) );
605 }
606 
607 
608 double QgsDistanceArea::getQbar( double x ) const
609 {
610  double cosx, cosx2;
611 
612  cosx = cos( x );
613  cosx2 = cosx * cosx;
614 
615  return cosx * ( m_QbarA + cosx2 * ( m_QbarB + cosx2 * ( m_QbarC + cosx2 * m_QbarD ) ) );
616 }
617 
618 
619 void QgsDistanceArea::computeAreaInit()
620 {
621  //don't try to perform calculations if no ellipsoid
622  if ( mEllipsoid == GEO_NONE )
623  {
624  return;
625  }
626 
627  double a2 = ( mSemiMajor * mSemiMajor );
628  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
629  double e4, e6;
630 
631  m_TwoPI = M_PI + M_PI;
632 
633  e4 = e2 * e2;
634  e6 = e4 * e2;
635 
636  m_AE = a2 * ( 1 - e2 );
637 
638  m_QA = ( 2.0 / 3.0 ) * e2;
639  m_QB = ( 3.0 / 5.0 ) * e4;
640  m_QC = ( 4.0 / 7.0 ) * e6;
641 
642  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
643  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
644  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
645  m_QbarD = ( 4.0 / 49.0 ) * e6;
646 
647  m_Qp = getQ( M_PI / 2 );
648  m_E = 4 * M_PI * m_Qp * m_AE;
649  if ( m_E < 0.0 )
650  m_E = -m_E;
651 }
652 
653 void QgsDistanceArea::setFromParams( const QgsEllipsoidUtils::EllipsoidParameters &params )
654 {
655  if ( params.useCustomParameters )
656  {
657  setEllipsoid( params.semiMajor, params.semiMinor );
658  }
659  else
660  {
661  mSemiMajor = params.semiMajor;
662  mSemiMinor = params.semiMinor;
663  mInvFlattening = params.inverseFlattening;
664  mCoordTransform.setDestinationCrs( params.crs );
665  // precalculate some values for area calculations
666  computeAreaInit();
667  }
668 }
669 
670 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint> &points ) const
671 {
672  if ( points.isEmpty() )
673  {
674  return 0;
675  }
676 
677  double x1, y1, x2, y2, dx, dy;
678  double Qbar1, Qbar2;
679  double area;
680 
681  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
682  if ( !willUseEllipsoid() )
683  {
684  return computePolygonFlatArea( points );
685  }
686  int n = points.size();
687  x2 = DEG2RAD( points[n - 1].x() );
688  y2 = DEG2RAD( points[n - 1].y() );
689  Qbar2 = getQbar( y2 );
690 
691  area = 0.0;
692 
693  for ( int i = 0; i < n; i++ )
694  {
695  x1 = x2;
696  y1 = y2;
697  Qbar1 = Qbar2;
698 
699  x2 = DEG2RAD( points[i].x() );
700  y2 = DEG2RAD( points[i].y() );
701  Qbar2 = getQbar( y2 );
702 
703  if ( x1 > x2 )
704  while ( x1 - x2 > M_PI )
705  x2 += m_TwoPI;
706  else if ( x2 > x1 )
707  while ( x2 - x1 > M_PI )
708  x1 += m_TwoPI;
709 
710  dx = x2 - x1;
711  area += dx * ( m_Qp - getQ( y2 ) );
712 
713  dy = y2 - y1;
714  if ( !qgsDoubleNear( dy, 0.0 ) )
715  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
716  }
717  if ( ( area *= m_AE ) < 0.0 )
718  area = -area;
719 
720  /* kludge - if polygon circles the south pole the area will be
721  * computed as if it cirlced the north pole. The correction is
722  * the difference between total surface area of the earth and
723  * the "north pole" area.
724  */
725  if ( area > m_E )
726  area = m_E;
727  if ( area > m_E / 2 )
728  area = m_E - area;
729 
730  return area;
731 }
732 
733 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint> &points ) const
734 {
735  // Normal plane area calculations.
736  double area = 0.0;
737  int i, size;
738 
739  size = points.size();
740 
741  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
742  for ( i = 0; i < size; i++ )
743  {
744  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
745  // Using '% size', so that we always end with the starting point
746  // and thus close the polygon.
747  area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
748  }
749  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
750  area = area / 2.0;
751  return qAbs( area ); // All areas are positive!
752 }
753 
754 QString QgsDistanceArea::formatDistance( double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit )
755 {
756  return QgsUnitTypes::formatDistance( distance, decimals, unit, keepBaseUnit );
757 }
758 
759 QString QgsDistanceArea::formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit )
760 {
761  return QgsUnitTypes::formatArea( area, decimals, unit, keepBaseUnit );
762 }
763 
765 {
766  // get the conversion factor between the specified units
767  QgsUnitTypes::DistanceUnit measureUnits = lengthUnits();
768  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
769 
770  double result = length * factorUnits;
771  QgsDebugMsgLevel( QString( "Converted length of %1 %2 to %3 %4" ).arg( length )
772  .arg( QgsUnitTypes::toString( measureUnits ) )
773  .arg( result )
774  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
775  return result;
776 }
777 
779 {
780  // get the conversion factor between the specified units
781  QgsUnitTypes::AreaUnit measureUnits = areaUnits();
782  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
783 
784  double result = area * factorUnits;
785  QgsDebugMsgLevel( QString( "Converted area of %1 %2 to %3 %4" ).arg( area )
786  .arg( QgsUnitTypes::toString( measureUnits ) )
787  .arg( result )
788  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
789  return result;
790 }
static void convertPointList(const QList< QgsPoint > &input, QgsPointSequence &output)
Upgrades a point list from QgsPoint to QgsPointV2.
const QgsCurve * interiorRing(int i) const
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
double y
Definition: qgspoint.h:42
QList< QgsPointV2 > QgsPointSequence
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.
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspoint.cpp:274
double measureLineProjected(const QgsPoint &p1, double distance=1, double azimuth=M_PI/2, QgsPoint *projectedPoint=nullptr) const
Calculates the distance from one point with distance in meters and azimuth (direction) When the sourc...
#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:193
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.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from...
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
QgsPoint transform(const QgsPoint &point, TransformDirection direction=ForwardTransform) const
Transform the point from the source CRS to the destination CRS.
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:181
QgsPoint computeSpheroidProject(const QgsPoint &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...
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
Computes the bearing (in radians) between two points.
#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:306
double measureLine(const QList< QgsPoint > &points) const
Measures the length of a line with multiple segments.
double measurePolygon(const QList< QgsPoint > &points) const
Measures the area of the polygon described by a set of points.
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
A class to represent a point.
Definition: qgspoint.h:37
#define POW2(x)
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:45
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
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.
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.
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.
void setSourceCrs(const QgsCoordinateReferenceSystem &srcCRS)
Sets source spatial reference system.
QgsPoint project(double distance, double bearing) const
Returns a new point which corresponds to this point projected by a specified distance in a specified ...
Definition: qgspoint.cpp:291
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.
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.
double x
Definition: qgspoint.h:41