QGIS API Documentation  3.14.0-Pi (9f7028fd23)
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 #include "qgsmultilinestring.h"
37 
38 #include <geodesic.h>
39 
40 #define DEG2RAD(x) ((x)*M_PI/180)
41 #define RAD2DEG(r) (180.0 * (r) / M_PI)
42 #define POW2(x) ((x)*(x))
43 
45 {
46  // init with default settings
47  mSemiMajor = -1.0;
48  mSemiMinor = -1.0;
49  mInvFlattening = -1.0;
50  QgsCoordinateTransformContext context; // this is ok - by default we have a source/dest of WGS84, so no reprojection takes place
51  setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), context ); // WGS 84
52  setEllipsoid( geoNone() );
53 }
54 
56 {
57  return mEllipsoid != geoNone();
58 }
59 
61 {
62  mCoordTransform.setContext( context );
63  mCoordTransform.setSourceCrs( srcCRS );
64 }
65 
66 bool QgsDistanceArea::setEllipsoid( const QString &ellipsoid )
67 {
68  // Shortcut if ellipsoid is none.
69  if ( ellipsoid == geoNone() )
70  {
71  mEllipsoid = geoNone();
72  return true;
73  }
74 
76  if ( !params.valid )
77  {
78  return false;
79  }
80  else
81  {
82  mEllipsoid = ellipsoid;
83  setFromParams( params );
84  return true;
85  }
86 }
87 
88 // Inverse flattening is calculated with invf = a/(a-b)
89 // Also, b = a-(a/invf)
90 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
91 {
92  mEllipsoid = QStringLiteral( "PARAMETER:%1:%2" ).arg( qgsDoubleToString( semiMajor ), qgsDoubleToString( semiMinor ) );
93  mSemiMajor = semiMajor;
94  mSemiMinor = semiMinor;
95  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
96 
97  computeAreaInit();
98 
99  return true;
100 }
101 
102 double QgsDistanceArea::measure( const QgsAbstractGeometry *geomV2, MeasureType type ) const
103 {
104  if ( !geomV2 )
105  {
106  return 0.0;
107  }
108 
109  int geomDimension = geomV2->dimension();
110  if ( geomDimension <= 0 )
111  {
112  return 0.0;
113  }
114 
115  MeasureType measureType = type;
116  if ( measureType == Default )
117  {
118  measureType = ( geomDimension == 1 ? Length : Area );
119  }
120 
121  if ( !willUseEllipsoid() )
122  {
123  //no transform required
124  if ( measureType == Length )
125  {
126  return geomV2->length();
127  }
128  else
129  {
130  return geomV2->area();
131  }
132  }
133  else
134  {
135  //multigeom is sum of measured parts
136  const QgsGeometryCollection *collection = qgsgeometry_cast<const QgsGeometryCollection *>( geomV2 );
137  if ( collection )
138  {
139  double sum = 0;
140  for ( int i = 0; i < collection->numGeometries(); ++i )
141  {
142  sum += measure( collection->geometryN( i ), measureType );
143  }
144  return sum;
145  }
146 
147  if ( measureType == Length )
148  {
149  const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( geomV2 );
150  if ( !curve )
151  {
152  return 0.0;
153  }
154 
155  QgsLineString *lineString = curve->curveToLine();
156  double length = measureLine( lineString );
157  delete lineString;
158  return length;
159  }
160  else
161  {
162  const QgsSurface *surface = qgsgeometry_cast<const QgsSurface *>( geomV2 );
163  if ( !surface )
164  return 0.0;
165 
166  QgsPolygon *polygon = surface->surfaceToPolygon();
167 
168  double area = 0;
169  const QgsCurve *outerRing = polygon->exteriorRing();
170  area += measurePolygon( outerRing );
171 
172  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
173  {
174  const QgsCurve *innerRing = polygon->interiorRing( i );
175  area -= measurePolygon( innerRing );
176  }
177  delete polygon;
178  return area;
179  }
180  }
181 }
182 
183 double QgsDistanceArea::measureArea( const QgsGeometry &geometry ) const
184 {
185  if ( geometry.isNull() )
186  return 0.0;
187 
188  const QgsAbstractGeometry *geomV2 = geometry.constGet();
189  return measure( geomV2, Area );
190 }
191 
192 double QgsDistanceArea::measureLength( const QgsGeometry &geometry ) const
193 {
194  if ( geometry.isNull() )
195  return 0.0;
196 
197  const QgsAbstractGeometry *geomV2 = geometry.constGet();
198  return measure( geomV2, Length );
199 }
200 
201 double QgsDistanceArea::measurePerimeter( const QgsGeometry &geometry ) const
202 {
203  if ( geometry.isNull() )
204  return 0.0;
205 
206  const QgsAbstractGeometry *geomV2 = geometry.constGet();
207  if ( !geomV2 || geomV2->dimension() < 2 )
208  {
209  return 0.0;
210  }
211 
212  if ( !willUseEllipsoid() )
213  {
214  return geomV2->perimeter();
215  }
216 
217  //create list with (single) surfaces
218  QVector< const QgsSurface * > surfaces;
219  const QgsSurface *surf = qgsgeometry_cast<const QgsSurface *>( geomV2 );
220  if ( surf )
221  {
222  surfaces.append( surf );
223  }
224  const QgsMultiSurface *multiSurf = qgsgeometry_cast<const QgsMultiSurface *>( geomV2 );
225  if ( multiSurf )
226  {
227  surfaces.reserve( ( surf ? 1 : 0 ) + multiSurf->numGeometries() );
228  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
229  {
230  surfaces.append( static_cast<const QgsSurface *>( multiSurf->geometryN( i ) ) );
231  }
232  }
233 
234  double length = 0;
235  QVector<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
236  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
237  {
238  if ( !*surfaceIt )
239  {
240  continue;
241  }
242 
243  QgsPolygon *poly = ( *surfaceIt )->surfaceToPolygon();
244  const QgsCurve *outerRing = poly->exteriorRing();
245  if ( outerRing )
246  {
247  length += measure( outerRing );
248  }
249  int nInnerRings = poly->numInteriorRings();
250  for ( int i = 0; i < nInnerRings; ++i )
251  {
252  length += measure( poly->interiorRing( i ) );
253  }
254  delete poly;
255  }
256  return length;
257 }
258 
259 double QgsDistanceArea::measureLine( const QgsCurve *curve ) const
260 {
261  if ( !curve )
262  {
263  return 0.0;
264  }
265 
266  QgsPointSequence linePointsV2;
267  QVector<QgsPointXY> linePoints;
268  curve->points( linePointsV2 );
269  QgsGeometry::convertPointList( linePointsV2, linePoints );
270  return measureLine( linePoints );
271 }
272 
273 double QgsDistanceArea::measureLine( const QVector<QgsPointXY> &points ) const
274 {
275  if ( points.size() < 2 )
276  return 0;
277 
278  double total = 0;
279  QgsPointXY p1, p2;
280 
281  try
282  {
283  if ( willUseEllipsoid() )
284  p1 = mCoordTransform.transform( points[0] );
285  else
286  p1 = points[0];
287 
288  for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
289  {
290  if ( willUseEllipsoid() )
291  {
292  p2 = mCoordTransform.transform( *i );
293  total += computeDistanceBearing( p1, p2 );
294  }
295  else
296  {
297  p2 = *i;
298  total += measureLine( p1, p2 );
299  }
300 
301  p1 = p2;
302  }
303 
304  return total;
305  }
306  catch ( QgsCsException &cse )
307  {
308  Q_UNUSED( cse )
309  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
310  return 0.0;
311  }
312 
313 }
314 
315 double QgsDistanceArea::measureLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const
316 {
317  double result;
318 
319  try
320  {
321  QgsPointXY pp1 = p1, pp2 = p2;
322 
323  QgsDebugMsgLevel( QStringLiteral( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
324  if ( willUseEllipsoid() )
325  {
326  QgsDebugMsgLevel( QStringLiteral( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
327  QgsDebugMsgLevel( QStringLiteral( "From proj4 : %1" ).arg( mCoordTransform.sourceCrs().toProj() ), 4 );
328  QgsDebugMsgLevel( QStringLiteral( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj() ), 4 );
329  pp1 = mCoordTransform.transform( p1 );
330  pp2 = mCoordTransform.transform( p2 );
331  QgsDebugMsgLevel( QStringLiteral( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
332  result = computeDistanceBearing( pp1, pp2 );
333  }
334  else
335  {
336  QgsDebugMsgLevel( QStringLiteral( "Cartesian calculation on canvas coordinates" ), 4 );
337  result = p2.distance( p1 );
338  }
339  }
340  catch ( QgsCsException &cse )
341  {
342  Q_UNUSED( cse )
343  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
344  result = 0.0;
345  }
346  QgsDebugMsgLevel( QStringLiteral( "The result was %1" ).arg( result ), 3 );
347  return result;
348 }
349 
350 double QgsDistanceArea::measureLineProjected( const QgsPointXY &p1, double distance, double azimuth, QgsPointXY *projectedPoint ) const
351 {
352  double result = 0.0;
353  QgsPointXY p2;
354  if ( mCoordTransform.sourceCrs().isGeographic() && willUseEllipsoid() )
355  {
356  p2 = computeSpheroidProject( p1, distance, azimuth );
357  result = p1.distance( p2 );
358  }
359  else // Cartesian coordinates
360  {
361  result = distance; // Avoid rounding errors when using meters [return as sent]
362  if ( sourceCrs().mapUnits() != QgsUnitTypes::DistanceMeters )
363  {
364  distance = ( distance * QgsUnitTypes::fromUnitToUnitFactor( QgsUnitTypes::DistanceMeters, sourceCrs().mapUnits() ) );
365  result = p1.distance( p2 );
366  }
367  p2 = p1.project( distance, azimuth );
368  }
369  QgsDebugMsgLevel( 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]" )
370  .arg( QString::number( distance, 'f', 7 ),
372  QString::number( result, 'f', 7 ),
373  mCoordTransform.sourceCrs().isGeographic() ? QStringLiteral( "Geographic" ) : QStringLiteral( "Cartesian" ),
374  QgsUnitTypes::toString( sourceCrs().mapUnits() ) )
375  .arg( azimuth )
376  .arg( p1.asWkt(),
377  p2.asWkt(),
378  sourceCrs().description(),
379  mEllipsoid )
380  .arg( sourceCrs().isGeographic() )
381  .arg( QStringLiteral( "SemiMajor[%1] SemiMinor[%2] InvFlattening[%3] " ).arg( QString::number( mSemiMajor, 'f', 7 ), QString::number( mSemiMinor, 'f', 7 ), QString::number( mInvFlattening, 'f', 7 ) ) ), 4 );
382  if ( projectedPoint )
383  {
384  *projectedPoint = QgsPointXY( p2 );
385  }
386  return result;
387 }
388 
389 /*
390  * From original rttopo documentation:
391  * Tested against:
392  * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
393  * and
394  * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
395  */
397  const QgsPointXY &p1, double distance, double azimuth ) const
398 {
399  // ellipsoid
400  double a = mSemiMajor;
401  double b = mSemiMinor;
402  double f = 1 / mInvFlattening;
403  if ( ( ( a < 0 ) && ( b < 0 ) ) ||
404  ( ( p1.x() < -180.0 ) || ( p1.x() > 180.0 ) || ( p1.y() < -85.05115 ) || ( p1.y() > 85.05115 ) ) )
405  {
406  // latitudes outside these bounds cause the calculations to become unstable and can return invalid results
407  return QgsPoint( 0, 0 );
408 
409  }
410  double radians_lat = DEG2RAD( p1.y() );
411  double radians_long = DEG2RAD( p1.x() );
412  double b2 = POW2( b ); // spheroid_mu2
413  double omf = 1 - f;
414  double tan_u1 = omf * std::tan( radians_lat );
415  double u1 = std::atan( tan_u1 );
416  double sigma, last_sigma, delta_sigma, two_sigma_m;
417  double sigma1, sin_alpha, alpha, cos_alphasq;
418  double u2, A, B;
419  double lat2, lambda, lambda2, C, omega;
420  int i = 0;
421  if ( azimuth < 0.0 )
422  {
423  azimuth = azimuth + M_PI * 2.0;
424  }
425  if ( azimuth > ( M_PI * 2.0 ) )
426  {
427  azimuth = azimuth - M_PI * 2.0;
428  }
429  sigma1 = std::atan2( tan_u1, std::cos( azimuth ) );
430  sin_alpha = std::cos( u1 ) * std::sin( azimuth );
431  alpha = std::asin( sin_alpha );
432  cos_alphasq = 1.0 - POW2( sin_alpha );
433  u2 = POW2( std::cos( alpha ) ) * ( POW2( a ) - b2 ) / b2; // spheroid_mu2
434  A = 1.0 + ( u2 / 16384.0 ) * ( 4096.0 + u2 * ( -768.0 + u2 * ( 320.0 - 175.0 * u2 ) ) );
435  B = ( u2 / 1024.0 ) * ( 256.0 + u2 * ( -128.0 + u2 * ( 74.0 - 47.0 * u2 ) ) );
436  sigma = ( distance / ( b * A ) );
437  do
438  {
439  two_sigma_m = 2.0 * sigma1 + sigma;
440  delta_sigma = B * 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 ) ) ) ) ) );
441  last_sigma = sigma;
442  sigma = ( distance / ( b * A ) ) + delta_sigma;
443  i++;
444  }
445  while ( i < 999 && std::fabs( ( last_sigma - sigma ) / sigma ) > 1.0e-9 );
446 
447  lat2 = std::atan2( ( std::sin( u1 ) * std::cos( sigma ) + std::cos( u1 ) * std::sin( sigma ) *
448  std::cos( azimuth ) ), ( omf * std::sqrt( POW2( sin_alpha ) +
449  POW2( std::sin( u1 ) * std::sin( sigma ) - std::cos( u1 ) * std::cos( sigma ) *
450  std::cos( azimuth ) ) ) ) );
451  lambda = std::atan2( ( std::sin( sigma ) * std::sin( azimuth ) ), ( std::cos( u1 ) * std::cos( sigma ) -
452  std::sin( u1 ) * std::sin( sigma ) * std::cos( azimuth ) ) );
453  C = ( f / 16.0 ) * cos_alphasq * ( 4.0 + f * ( 4.0 - 3.0 * cos_alphasq ) );
454  omega = lambda - ( 1.0 - C ) * f * sin_alpha * ( sigma + C * std::sin( sigma ) *
455  ( std::cos( two_sigma_m ) + C * std::cos( sigma ) * ( -1.0 + 2.0 * POW2( std::cos( two_sigma_m ) ) ) ) );
456  lambda2 = radians_long + omega;
457  return QgsPointXY( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
458 }
459 
460 double QgsDistanceArea::latitudeGeodesicCrossesAntimeridian( const QgsPointXY &pp1, const QgsPointXY &pp2, double &fractionAlongLine ) const
461 {
462  QgsPointXY p1 = pp1;
463  QgsPointXY p2 = pp2;
464  if ( p1.x() < -120 )
465  p1.setX( p1.x() + 360 );
466  if ( p2.x() < -120 )
467  p2.setX( p2.x() + 360 );
468 
469  // we need p2.x() > 180 and p1.x() < 180
470  double p1x = p1.x() < 180 ? p1.x() : p2.x();
471  double p1y = p1.x() < 180 ? p1.y() : p2.y();
472  double p2x = p1.x() < 180 ? p2.x() : p1.x();
473  double p2y = p1.x() < 180 ? p2.y() : p1.y();
474  // lat/lon are our candidate intersection position - we want this to get as close to 180 as possible
475  // the first candidate is p2
476  double lat = p2y;
477  double lon = p2x;
478 
479  if ( mEllipsoid == geoNone() )
480  {
481  fractionAlongLine = ( 180 - p1x ) / ( p2x - p1x );
482  if ( p1.x() >= 180 )
483  fractionAlongLine = 1 - fractionAlongLine;
484  return p1y + ( 180 - p1x ) / ( p2x - p1x ) * ( p2y - p1y );
485  }
486 
487  geod_geodesic geod;
488  geod_init( &geod, mSemiMajor, 1 / mInvFlattening );
489 
490  geod_geodesicline line;
491  geod_inverseline( &line, &geod, p1y, p1x, p2y, p2x, GEOD_ALL );
492 
493  const double totalDist = line.s13;
494  double intersectionDist = line.s13;
495 
496  int iterations = 0;
497  double t = 0;
498  // iterate until our intersection candidate is within ~1 mm of the antimeridian (or too many iterations happened)
499  while ( std::fabs( lon - 180.0 ) > 0.00000001 && iterations < 100 )
500  {
501  if ( iterations > 0 && std::fabs( p2x - p1x ) > 5 )
502  {
503  // if we have too large a range of longitudes, we use a binary search to narrow the window -- this ensures we will converge
504  if ( lon < 180 )
505  {
506  p1x = lon;
507  p1y = lat;
508  }
509  else
510  {
511  p2x = lon;
512  p2y = lat;
513  }
514  QgsDebugMsgLevel( QStringLiteral( "Narrowed window to %1, %2 - %3, %4" ).arg( p1x ).arg( p1y ).arg( p2x ).arg( p2y ), 4 );
515 
516  geod_inverseline( &line, &geod, p1y, p1x, p2y, p2x, GEOD_ALL );
517  intersectionDist = line.s13 * 0.5;
518  }
519  else
520  {
521  // we have a sufficiently narrow window -- use Newton's method
522  // adjust intersection distance by fraction of how close the previous candidate was to 180 degrees longitude -
523  // this helps us close in to the correct longitude quickly
524  intersectionDist *= ( 180.0 - p1x ) / ( lon - p1x );
525  }
526 
527  // now work out the point on the geodesic this far from p1 - this becomes our new candidate for crossing the antimeridian
528 
529  geod_position( &line, intersectionDist, &lat, &lon, &t );
530  // we don't want to wrap longitudes > 180 around)
531  if ( lon < 0 )
532  lon += 360;
533 
534  iterations++;
535  QgsDebugMsgLevel( QStringLiteral( "After %1 iterations lon is %2, lat is %3, dist from p1: %4" ).arg( iterations ).arg( lon ).arg( lat ).arg( intersectionDist ), 4 );
536  }
537 
538  fractionAlongLine = intersectionDist / totalDist;
539  if ( p1.x() >= 180 )
540  fractionAlongLine = 1 - fractionAlongLine;
541 
542  // either converged on 180 longitude or hit too many iterations
543  return lat;
544 }
545 
547 {
549  return geometry;
550 
551  QgsGeometry g = geometry;
552  // TODO - avoid segmentization of curved geometries (if this is even possible!)
553  if ( QgsWkbTypes::isCurvedType( g.wkbType() ) )
555 
556  std::unique_ptr< QgsMultiLineString > res = qgis::make_unique< QgsMultiLineString >();
557  for ( auto part = g.const_parts_begin(); part != g.const_parts_end(); ++part )
558  {
559  const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( *part );
560  if ( !line )
561  continue;
562  if ( line->isEmpty() )
563  {
564  continue;
565  }
566 
567  std::unique_ptr< QgsLineString > l = qgis::make_unique< QgsLineString >();
568  try
569  {
570  double x = 0;
571  double y = 0;
572  double z = 0;
573  double m = 0;
574  QVector< QgsPoint > newPoints;
575  newPoints.reserve( line->numPoints() );
576  double prevLon = 0;
577  double prevLat = 0;
578  double lon = 0;
579  double lat = 0;
580  double prevZ = 0;
581  double prevM = 0;
582  for ( int i = 0; i < line->numPoints(); i++ )
583  {
584  QgsPoint p = line->pointN( i );
585  x = p.x();
586  if ( mCoordTransform.sourceCrs().isGeographic() )
587  {
588  x = std::fmod( x, 360.0 );
589  if ( x > 180 )
590  x -= 360;
591  p.setX( x );
592  }
593  y = p.y();
594  lon = x;
595  lat = y;
596  mCoordTransform.transformInPlace( lon, lat, z );
597 
598  //test if we crossed the antimeridian in this segment
599  if ( i > 0 && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
600  {
601  // we did!
602  // when crossing the antimeridian, we need to calculate the latitude
603  // at which the geodesic intersects the antimeridian
604  double fract = 0;
605  double lat180 = latitudeGeodesicCrossesAntimeridian( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ), fract );
606  if ( line->is3D() )
607  {
608  z = prevZ + ( p.z() - prevZ ) * fract;
609  }
610  if ( line->isMeasure() )
611  {
612  m = prevM + ( p.m() - prevM ) * fract;
613  }
614 
615  QgsPointXY antiMeridianPoint;
616  if ( prevLon < -120 )
617  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
618  else
619  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
620 
621  QgsPoint newPoint( antiMeridianPoint );
622  if ( line->is3D() )
623  newPoint.addZValue( z );
624  if ( line->isMeasure() )
625  newPoint.addMValue( m );
626 
627  if ( std::isfinite( newPoint.x() ) && std::isfinite( newPoint.y() ) )
628  {
629  newPoints << newPoint;
630  }
631  res->addGeometry( new QgsLineString( newPoints ) );
632 
633  newPoints.clear();
634  newPoints.reserve( line->numPoints() - i + 1 );
635 
636  if ( lon < -120 )
637  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
638  else
639  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
640 
641  if ( std::isfinite( antiMeridianPoint.x() ) && std::isfinite( antiMeridianPoint.y() ) )
642  {
643  // we want to keep the previously calculated z/m value for newPoint, if present. They're the same each
644  // of the antimeridian split
645  newPoint.setX( antiMeridianPoint.x() );
646  newPoint.setY( antiMeridianPoint.y() );
647  newPoints << newPoint;
648  }
649  }
650  newPoints << p;
651 
652  prevLon = lon;
653  prevLat = lat;
654  if ( line->is3D() )
655  prevZ = p.z();
656  if ( line->isMeasure() )
657  prevM = p.m();
658  }
659  res->addGeometry( new QgsLineString( newPoints ) );
660  }
661  catch ( QgsCsException & )
662  {
663  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform linestring. Unable to calculate break point." ) );
664  res->addGeometry( line->clone() );
665  break;
666  }
667  }
668 
669  return QgsGeometry( std::move( res ) );
670 }
671 
672 
673 QVector< QVector<QgsPointXY> > QgsDistanceArea::geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, const double interval, const bool breakLine ) const
674 {
675  if ( !willUseEllipsoid() )
676  {
677  return QVector< QVector< QgsPointXY > >() << ( QVector< QgsPointXY >() << p1 << p2 );
678  }
679 
680  geod_geodesic geod;
681  geod_init( &geod, mSemiMajor, 1 / mInvFlattening );
682 
683  QgsPointXY pp1, pp2;
684  try
685  {
686  pp1 = mCoordTransform.transform( p1 );
687  pp2 = mCoordTransform.transform( p2 );
688  }
689  catch ( QgsCsException & )
690  {
691  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate geodesic line." ) );
692  return QVector< QVector< QgsPointXY > >();
693  }
694 
695  geod_geodesicline line;
696  geod_inverseline( &line, &geod, pp1.y(), pp1.x(), pp2.y(), pp2.x(), GEOD_ALL );
697  const double totalDist = line.s13;
698 
699  QVector< QVector< QgsPointXY > > res;
700  QVector< QgsPointXY > currentPart;
701  currentPart << p1;
702  double d = interval;
703  double prevLon = pp1.x();
704  double prevLat = pp1.y();
705  bool lastRun = false;
706  double t = 0;
707  while ( true )
708  {
709  double lat, lon;
710  if ( lastRun )
711  {
712  lat = pp2.y();
713  lon = pp2.x();
714  if ( lon > 180 )
715  lon -= 360;
716  }
717  else
718  {
719  geod_position( &line, d, &lat, &lon, &t );
720  }
721 
722  if ( breakLine && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
723  {
724  // when breaking the geodesic at the antimeridian, we need to calculate the latitude
725  // at which the geodesic intersects the antimeridian, and add points to both line segments at this latitude
726  // on the antimeridian.
727  double fraction;
728  double lat180 = latitudeGeodesicCrossesAntimeridian( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ), fraction );
729 
730  try
731  {
732  QgsPointXY p;
733  if ( prevLon < -120 )
734  p = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
735  else
736  p = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
737 
738  if ( std::isfinite( p.x() ) && std::isfinite( p.y() ) )
739  currentPart << p;
740  }
741  catch ( QgsCsException & )
742  {
743  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point." ) );
744  }
745 
746  res << currentPart;
747  currentPart.clear();
748  try
749  {
750  QgsPointXY p;
751  if ( lon < -120 )
752  p = mCoordTransform.transform( QgsPointXY( -180, lat180 ), QgsCoordinateTransform::ReverseTransform );
753  else
754  p = mCoordTransform.transform( QgsPointXY( 180, lat180 ), QgsCoordinateTransform::ReverseTransform );
755 
756  if ( std::isfinite( p.x() ) && std::isfinite( p.y() ) )
757  currentPart << p;
758  }
759  catch ( QgsCsException & )
760  {
761  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point." ) );
762  }
763 
764  }
765 
766  prevLon = lon;
767  prevLat = lat;
768 
769  try
770  {
771  currentPart << mCoordTransform.transform( QgsPointXY( lon, lat ), QgsCoordinateTransform::ReverseTransform );
772  }
773  catch ( QgsCsException & )
774  {
775  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point." ) );
776  }
777 
778  if ( lastRun )
779  break;
780 
781  d += interval;
782  if ( d >= totalDist )
783  lastRun = true;
784  }
785  res << currentPart;
786  return res;
787 }
788 
790 {
791  return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();
792 }
793 
795 {
797  QgsUnitTypes::distanceToAreaUnit( mCoordTransform.sourceCrs().mapUnits() );
798 }
799 
800 double QgsDistanceArea::measurePolygon( const QgsCurve *curve ) const
801 {
802  if ( !curve )
803  {
804  return 0.0;
805  }
806 
807  QgsPointSequence linePointsV2;
808  curve->points( linePointsV2 );
809  QVector<QgsPointXY> linePoints;
810  QgsGeometry::convertPointList( linePointsV2, linePoints );
811  return measurePolygon( linePoints );
812 }
813 
814 
815 double QgsDistanceArea::measurePolygon( const QVector<QgsPointXY> &points ) const
816 {
817  try
818  {
819  if ( willUseEllipsoid() )
820  {
821  QVector<QgsPointXY> pts;
822  for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
823  {
824  pts.append( mCoordTransform.transform( *i ) );
825  }
826  return computePolygonArea( pts );
827  }
828  else
829  {
830  return computePolygonArea( points );
831  }
832  }
833  catch ( QgsCsException &cse )
834  {
835  Q_UNUSED( cse )
836  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
837  return 0.0;
838  }
839 }
840 
841 
842 double QgsDistanceArea::bearing( const QgsPointXY &p1, const QgsPointXY &p2 ) const
843 {
844  QgsPointXY pp1 = p1, pp2 = p2;
845  double bearing;
846 
847  if ( willUseEllipsoid() )
848  {
849  pp1 = mCoordTransform.transform( p1 );
850  pp2 = mCoordTransform.transform( p2 );
851  computeDistanceBearing( pp1, pp2, &bearing );
852  }
853  else //compute simple planar azimuth
854  {
855  double dx = p2.x() - p1.x();
856  double dy = p2.y() - p1.y();
857  bearing = std::atan2( dx, dy );
858  }
859 
860  return bearing;
861 }
862 
863 
865 // distance calculation
866 
867 double QgsDistanceArea::computeDistanceBearing(
868  const QgsPointXY &p1, const QgsPointXY &p2,
869  double *course1, double *course2 ) const
870 {
871  if ( qgsDoubleNear( p1.x(), p2.x() ) && qgsDoubleNear( p1.y(), p2.y() ) )
872  return 0;
873 
874  // ellipsoid
875  double a = mSemiMajor;
876  double b = mSemiMinor;
877  double f = 1 / mInvFlattening;
878 
879  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
880  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
881 
882  double L = p2_lon - p1_lon;
883  double U1 = std::atan( ( 1 - f ) * std::tan( p1_lat ) );
884  double U2 = std::atan( ( 1 - f ) * std::tan( p2_lat ) );
885  double sinU1 = std::sin( U1 ), cosU1 = std::cos( U1 );
886  double sinU2 = std::sin( U2 ), cosU2 = std::cos( U2 );
887  double lambda = L;
888  double lambdaP = 2 * M_PI;
889 
890  double sinLambda = 0;
891  double cosLambda = 0;
892  double sinSigma = 0;
893  double cosSigma = 0;
894  double sigma = 0;
895  double alpha = 0;
896  double cosSqAlpha = 0;
897  double cos2SigmaM = 0;
898  double C = 0;
899  double tu1 = 0;
900  double tu2 = 0;
901 
902  int iterLimit = 20;
903  while ( std::fabs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
904  {
905  sinLambda = std::sin( lambda );
906  cosLambda = std::cos( lambda );
907  tu1 = ( cosU2 * sinLambda );
908  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
909  sinSigma = std::sqrt( tu1 * tu1 + tu2 * tu2 );
910  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
911  sigma = std::atan2( sinSigma, cosSigma );
912  alpha = std::asin( cosU1 * cosU2 * sinLambda / sinSigma );
913  cosSqAlpha = std::cos( alpha ) * std::cos( alpha );
914  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
915  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
916  lambdaP = lambda;
917  lambda = L + ( 1 - C ) * f * std::sin( alpha ) *
918  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
919  }
920 
921  if ( iterLimit == 0 )
922  return -1; // formula failed to converge
923 
924  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
925  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
926  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
927  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
928  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
929  double s = b * A * ( sigma - deltaSigma );
930 
931  if ( course1 )
932  {
933  *course1 = std::atan2( tu1, tu2 );
934  }
935  if ( course2 )
936  {
937  // PI is added to return azimuth from P2 to P1
938  *course2 = std::atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
939  }
940 
941  return s;
942 }
943 
945 // stuff for measuring areas - copied from GRASS
946 // don't know how does it work, but it's working .)
947 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
948 
949 double QgsDistanceArea::getQ( double x ) const
950 {
951  double sinx, sinx2;
952 
953  sinx = std::sin( x );
954  sinx2 = sinx * sinx;
955 
956  return sinx * ( 1 + sinx2 * ( m_QA + sinx2 * ( m_QB + sinx2 * m_QC ) ) );
957 }
958 
959 
960 double QgsDistanceArea::getQbar( double x ) const
961 {
962  double cosx, cosx2;
963 
964  cosx = std::cos( x );
965  cosx2 = cosx * cosx;
966 
967  return cosx * ( m_QbarA + cosx2 * ( m_QbarB + cosx2 * ( m_QbarC + cosx2 * m_QbarD ) ) );
968 }
969 
970 
971 void QgsDistanceArea::computeAreaInit()
972 {
973  //don't try to perform calculations if no ellipsoid
974  if ( mEllipsoid == geoNone() )
975  {
976  return;
977  }
978 
979  double a2 = ( mSemiMajor * mSemiMajor );
980  double e2 = 1 - ( ( mSemiMinor * mSemiMinor ) / a2 );
981  double e4, e6;
982 
983  m_TwoPI = M_PI + M_PI;
984 
985  e4 = e2 * e2;
986  e6 = e4 * e2;
987 
988  m_AE = a2 * ( 1 - e2 );
989 
990  m_QA = ( 2.0 / 3.0 ) * e2;
991  m_QB = ( 3.0 / 5.0 ) * e4;
992  m_QC = ( 4.0 / 7.0 ) * e6;
993 
994  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
995  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
996  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
997  m_QbarD = ( 4.0 / 49.0 ) * e6;
998 
999  m_Qp = getQ( M_PI_2 );
1000  m_E = 4 * M_PI * m_Qp * m_AE;
1001  if ( m_E < 0.0 )
1002  m_E = -m_E;
1003 }
1004 
1005 void QgsDistanceArea::setFromParams( const QgsEllipsoidUtils::EllipsoidParameters &params )
1006 {
1007  if ( params.useCustomParameters )
1008  {
1009  setEllipsoid( params.semiMajor, params.semiMinor );
1010  }
1011  else
1012  {
1013  mSemiMajor = params.semiMajor;
1014  mSemiMinor = params.semiMinor;
1015  mInvFlattening = params.inverseFlattening;
1016  mCoordTransform.setDestinationCrs( params.crs );
1017  // precalculate some values for area calculations
1018  computeAreaInit();
1019  }
1020 }
1021 
1022 double QgsDistanceArea::computePolygonArea( const QVector<QgsPointXY> &points ) const
1023 {
1024  if ( points.isEmpty() )
1025  {
1026  return 0;
1027  }
1028 
1029  // IMPORTANT
1030  // don't change anything here without reporting the changes to upstream (GRASS)
1031  // let's all be good opensource citizens and share the improvements!
1032 
1033  double x1, y1, x2, y2, dx, dy;
1034  double Qbar1, Qbar2;
1035  double area;
1036 
1037  /* GRASS comment: threshold for dy, should be between 1e-4 and 1e-7
1038  * See relevant discussion at https://trac.osgeo.org/grass/ticket/3369
1039  */
1040  const double thresh = 1e-6;
1041 
1042  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
1043  if ( !willUseEllipsoid() )
1044  {
1045  return computePolygonFlatArea( points );
1046  }
1047  int n = points.size();
1048  x2 = DEG2RAD( points[n - 1].x() );
1049  y2 = DEG2RAD( points[n - 1].y() );
1050  Qbar2 = getQbar( y2 );
1051 
1052  area = 0.0;
1053 
1054  for ( int i = 0; i < n; i++ )
1055  {
1056  x1 = x2;
1057  y1 = y2;
1058  Qbar1 = Qbar2;
1059 
1060  x2 = DEG2RAD( points[i].x() );
1061  y2 = DEG2RAD( points[i].y() );
1062  Qbar2 = getQbar( y2 );
1063 
1064  if ( x1 > x2 )
1065  while ( x1 - x2 > M_PI )
1066  x2 += m_TwoPI;
1067  else if ( x2 > x1 )
1068  while ( x2 - x1 > M_PI )
1069  x1 += m_TwoPI;
1070 
1071  dx = x2 - x1;
1072  dy = y2 - y1;
1073  if ( std::fabs( dy ) > thresh )
1074  {
1075  /* account for different latitudes y1, y2 */
1076  area += dx * ( m_Qp - ( Qbar2 - Qbar1 ) / dy );
1077  }
1078  else
1079  {
1080  /* latitudes y1, y2 are (nearly) identical */
1081 
1082  /* if y2 becomes similar to y1, i.e. y2 -> y1
1083  * Qbar2 - Qbar1 -> 0 and dy -> 0
1084  * (Qbar2 - Qbar1) / dy -> ?
1085  * (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
1086  * Metz 2017
1087  */
1088  area += dx * ( m_Qp - getQ( ( y1 + y2 ) / 2.0 ) );
1089 
1090  /* original:
1091  * area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
1092  */
1093  }
1094  }
1095  if ( ( area *= m_AE ) < 0.0 )
1096  area = -area;
1097 
1098  /* kludge - if polygon circles the south pole the area will be
1099  * computed as if it cirlced the north pole. The correction is
1100  * the difference between total surface area of the earth and
1101  * the "north pole" area.
1102  */
1103  if ( area > m_E )
1104  area = m_E;
1105  if ( area > m_E / 2 )
1106  area = m_E - area;
1107 
1108  return area;
1109 }
1110 
1111 double QgsDistanceArea::computePolygonFlatArea( const QVector<QgsPointXY> &points ) const
1112 {
1113  // Normal plane area calculations.
1114  double area = 0.0;
1115  int i, size;
1116 
1117  size = points.size();
1118 
1119  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
1120  for ( i = 0; i < size; i++ )
1121  {
1122  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
1123  // Using '% size', so that we always end with the starting point
1124  // and thus close the polygon.
1125  area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
1126  }
1127  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
1128  area = area / 2.0;
1129  return std::fabs( area ); // All areas are positive!
1130 }
1131 
1132 QString QgsDistanceArea::formatDistance( double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit )
1133 {
1134  return QgsUnitTypes::formatDistance( distance, decimals, unit, keepBaseUnit );
1135 }
1136 
1137 QString QgsDistanceArea::formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit )
1138 {
1139  return QgsUnitTypes::formatArea( area, decimals, unit, keepBaseUnit );
1140 }
1141 
1143 {
1144  // get the conversion factor between the specified units
1145  QgsUnitTypes::DistanceUnit measureUnits = lengthUnits();
1146  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1147 
1148  double result = length * factorUnits;
1149  QgsDebugMsgLevel( QStringLiteral( "Converted length of %1 %2 to %3 %4" ).arg( length )
1150  .arg( QgsUnitTypes::toString( measureUnits ) )
1151  .arg( result )
1152  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
1153  return result;
1154 }
1155 
1157 {
1158  // get the conversion factor between the specified units
1159  QgsUnitTypes::AreaUnit measureUnits = areaUnits();
1160  double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
1161 
1162  double result = area * factorUnits;
1163  QgsDebugMsgLevel( QStringLiteral( "Converted area of %1 %2 to %3 %4" ).arg( area )
1164  .arg( QgsUnitTypes::toString( measureUnits ) )
1165  .arg( result )
1166  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
1167  return result;
1168 }
QgsCurve
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
POW2
#define POW2(x)
Definition: qgsdistancearea.cpp:42
qgspolygon.h
QgsDistanceArea::ellipsoid
QString ellipsoid() const
Returns ellipsoid's acronym.
Definition: qgsdistancearea.h:112
QgsPoint::setY
void setY(double y)
Sets the point's y-coordinate.
Definition: qgspoint.h:298
QgsLineString::pointN
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
Definition: qgslinestring.cpp:660
QgsDistanceArea::QgsDistanceArea
QgsDistanceArea()
Constructor.
Definition: qgsdistancearea.cpp:44
QgsPointXY::y
double y
Definition: qgspointxy.h:48
QgsDistanceArea::measureLength
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
Definition: qgsdistancearea.cpp:192
QgsCoordinateTransformContext
Definition: qgscoordinatetransformcontext.h:57
QgsPolygon
Polygon geometry type.
Definition: qgspolygon.h:33
qgslinestring.h
QgsPointXY::project
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:87
QgsPoint
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
QgsGeometry::const_parts_end
QgsAbstractGeometry::const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary part after the last part of the geometry.
Definition: qgsgeometry.cpp:1895
qgswkbptr.h
QgsEllipsoidUtils::EllipsoidParameters::useCustomParameters
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)
Definition: qgsellipsoidutils.h:50
QgsGeometryCollection::numGeometries
int numGeometries() const
Returns the number of geometries within the collection.
Definition: qgsgeometrycollection.h:51
QgsDebugMsgLevel
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
QgsPoint::addZValue
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:511
QgsLineString::numPoints
int numPoints() const override
Returns the number of points in the curve.
Definition: qgslinestring.cpp:650
QgsCoordinateTransform::setContext
void setContext(const QgsCoordinateTransformContext &context)
Sets the context in which the coordinate transform should be calculated.
Definition: qgscoordinatetransform.cpp:200
QgsPoint::clear
void clear() override
Clears the geometry, ie reset it to a null geometry.
Definition: qgspoint.cpp:324
QgsUnitTypes::distanceToAreaUnit
static Q_INVOKABLE QgsUnitTypes::AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
Definition: qgsunittypes.cpp:1176
QgsCoordinateTransform::setSourceCrs
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
Definition: qgscoordinatetransform.cpp:159
QgsCurvePolygon::exteriorRing
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
Definition: qgscurvepolygon.h:86
QgsGeometry::wkbType
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Definition: qgsgeometry.cpp:345
QgsPoint::z
double z
Definition: qgspoint.h:60
QgsDistanceArea::sourceCrs
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
Definition: qgsdistancearea.h:76
qgis.h
QgsDistanceArea::measureArea
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
Definition: qgsdistancearea.cpp:183
QgsSurface
Definition: qgssurface.h:32
QgsAbstractGeometry::length
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
Definition: qgsabstractgeometry.cpp:132
qgsunittypes.h
QgsPointXY::asWkt
QString asWkt() const
Returns the well known text representation for the point (e.g.
Definition: qgspointxy.cpp:69
QgsGeometry::const_parts_begin
QgsAbstractGeometry::const_part_iterator const_parts_begin() const
Returns STL-style const iterator pointing to the first part of the geometry.
Definition: qgsgeometry.cpp:1888
QgsMultiSurface
Multi surface geometry collection.
Definition: qgsmultisurface.h:29
QgsEllipsoidUtils::ellipsoidParameters
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
Definition: qgsellipsoidutils.cpp:41
QgsDistanceArea::splitGeometryAtAntimeridian
QgsGeometry splitGeometryAtAntimeridian(const QgsGeometry &geometry) const
Splits a (Multi)LineString geometry at the antimeridian (longitude +/- 180 degrees).
Definition: qgsdistancearea.cpp:546
QgsUnitTypes::formatArea
static Q_INVOKABLE QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
Definition: qgsunittypes.cpp:2864
QgsLineString
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
QgsUnitTypes::DistanceUnit
DistanceUnit
Units of distance.
Definition: qgsunittypes.h:67
QgsCoordinateTransform::transform
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
Definition: qgscoordinatetransform.cpp:239
QgsCoordinateTransform::ReverseTransform
@ ReverseTransform
Transform from destination to source CRS.
Definition: qgscoordinatetransform.h:61
QgsCurvePolygon::interiorRing
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
Definition: qgscurvepolygon.h:99
QgsDistanceArea::measureLineProjected
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...
Definition: qgsdistancearea.cpp:350
qgsDoubleToString
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:275
QgsEllipsoidUtils::EllipsoidParameters::crs
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
Definition: qgsellipsoidutils.h:56
QgsDistanceArea::convertAreaMeasurement
double convertAreaMeasurement(double area, QgsUnitTypes::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
Definition: qgsdistancearea.cpp:1156
QgsDistanceArea::measurePolygon
double measurePolygon(const QVector< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
Definition: qgsdistancearea.cpp:815
QgsCoordinateTransform::destinationCrs
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
Definition: qgscoordinatetransform.cpp:234
QgsGeometry::convertToStraightSegment
void convertToStraightSegment(double tolerance=M_PI/180., QgsAbstractGeometry::SegmentationToleranceType toleranceType=QgsAbstractGeometry::MaximumAngle)
Converts the geometry to straight line segments, if it is a curved geometry type.
Definition: qgsgeometry.cpp:2815
QgsDistanceArea::areaUnits
QgsUnitTypes::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
Definition: qgsdistancearea.cpp:794
QgsUnitTypes::fromUnitToUnitFactor
static Q_INVOKABLE double fromUnitToUnitFactor(QgsUnitTypes::DistanceUnit fromUnit, QgsUnitTypes::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
Definition: qgsunittypes.cpp:352
QgsAbstractGeometry::area
virtual double area() const
Returns the planar, 2-dimensional area of the geometry.
Definition: qgsabstractgeometry.cpp:142
QgsUnitTypes::AreaSquareMeters
@ AreaSquareMeters
Square meters.
Definition: qgsunittypes.h:95
QgsSurface::surfaceToPolygon
virtual QgsPolygon * surfaceToPolygon() const =0
Gets a polygon representation of this surface.
QgsPoint::y
double y
Definition: qgspoint.h:59
QgsCoordinateTransform::sourceCrs
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from.
Definition: qgscoordinatetransform.cpp:229
QgsDistanceArea::setEllipsoid
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
Definition: qgsdistancearea.cpp:66
QgsGeometryCollection
Geometry collection.
Definition: qgsgeometrycollection.h:35
QgsCoordinateReferenceSystem::isGeographic
bool isGeographic
Definition: qgscoordinatereferencesystem.h:211
QgsDistanceArea::geodesicLine
QVector< QVector< QgsPointXY > > geodesicLine(const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine=false) const
Calculates the geodesic line between p1 and p2, which represents the shortest path on the ellipsoid b...
Definition: qgsdistancearea.cpp:673
QgsAbstractGeometry::dimension
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
QgsCsException
Definition: qgsexception.h:65
QgsLineString::clone
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
Definition: qgslinestring.cpp:293
QgsPoint::addMValue
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:522
QgsDistanceArea::bearing
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
Definition: qgsdistancearea.cpp:842
QgsCoordinateTransform::setDestinationCrs
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
Definition: qgscoordinatetransform.cpp:179
geoNone
CONSTLATIN1STRING geoNone()
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.h:672
QgsEllipsoidUtils::EllipsoidParameters::valid
bool valid
Whether ellipsoid parameters are valid.
Definition: qgsellipsoidutils.h:42
QgsUnitTypes::toString
static Q_INVOKABLE QString toString(QgsUnitTypes::DistanceUnit unit)
Returns a translated string representing a distance unit.
Definition: qgsunittypes.cpp:199
QgsPointXY::toString
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
Definition: qgspointxy.cpp:51
qgsDoubleNear
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:315
QgsPoint::setX
void setX(double x)
Sets the point's x-coordinate.
Definition: qgspoint.h:287
QgsPointXY::distance
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:196
QgsDistanceArea::measureLine
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
Definition: qgsdistancearea.cpp:273
QgsDistanceArea::lengthUnits
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
Definition: qgsdistancearea.cpp:789
QgsUnitTypes::formatDistance
static Q_INVOKABLE QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
Definition: qgsunittypes.cpp:2852
QgsWkbTypes::geometryType
static GeometryType geometryType(Type type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:937
QgsUnitTypes::DistanceMeters
@ DistanceMeters
Meters.
Definition: qgsunittypes.h:69
QgsPoint::m
double m
Definition: qgspoint.h:61
QgsDistanceArea::formatDistance
static QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
Definition: qgsdistancearea.cpp:1132
qgscoordinatetransform.h
QgsUnitTypes::AreaUnit
AreaUnit
Units of area.
Definition: qgsunittypes.h:93
QgsAbstractGeometry::is3D
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
Definition: qgsabstractgeometry.h:202
QgsGeometry::isNull
bool isNull
Definition: qgsgeometry.h:125
QgsGeometry::constGet
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Definition: qgsgeometry.cpp:128
QgsDistanceArea::setSourceCrs
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
Definition: qgsdistancearea.cpp:60
QgsEllipsoidUtils::EllipsoidParameters::semiMinor
double semiMinor
Semi-minor axis.
Definition: qgsellipsoidutils.h:47
QgsCoordinateReferenceSystem::toProj
QString toProj() const
Returns a Proj string representation of this CRS.
Definition: qgscoordinatereferencesystem.cpp:1398
QgsCoordinateReferenceSystem
Definition: qgscoordinatereferencesystem.h:206
QgsAbstractGeometry
Abstract base class for all geometries.
Definition: qgsabstractgeometry.h:71
DEG2RAD
#define DEG2RAD(x)
Definition: qgsdistancearea.cpp:40
QgsPointXY
Definition: qgspointxy.h:43
QgsLineString::isEmpty
bool isEmpty() const override
Returns true if the geometry is empty.
Definition: qgslinestring.cpp:308
QgsEllipsoidUtils::EllipsoidParameters
Contains parameters for an ellipsoid.
Definition: qgsellipsoidutils.h:39
QgsWkbTypes::LineGeometry
@ LineGeometry
Definition: qgswkbtypes.h:142
QgsCurvePolygon::numInteriorRings
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
Definition: qgscurvepolygon.h:76
QgsCoordinateReferenceSystem::mapUnits
QgsUnitTypes::DistanceUnit mapUnits
Definition: qgscoordinatereferencesystem.h:210
QgsGeometryCollection::geometryN
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
Definition: qgsgeometrycollection.h:79
qgsgeometry.h
QgsGeometry::convertPointList
static void convertPointList(const QVector< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPoint.
Definition: qgsgeometry.cpp:3043
QgsDistanceArea::measurePerimeter
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
Definition: qgsdistancearea.cpp:201
QgsPointSequence
QVector< QgsPoint > QgsPointSequence
Definition: qgsabstractgeometry.h:44
QgsGeometry
Definition: qgsgeometry.h:122
QgsDistanceArea::convertLengthMeasurement
double convertLengthMeasurement(double length, QgsUnitTypes::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
Definition: qgsdistancearea.cpp:1142
QgsPointXY::x
double x
Definition: qgspointxy.h:47
QgsPolygon::surfaceToPolygon
QgsPolygon * surfaceToPolygon() const override
Gets a polygon representation of this surface.
Definition: qgspolygon.cpp:306
RAD2DEG
#define RAD2DEG(r)
Definition: qgsdistancearea.cpp:41
QgsMessageLog::logMessage
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).
Definition: qgsmessagelog.cpp:27
QgsAbstractGeometry::perimeter
virtual double perimeter() const
Returns the planar, 2-dimensional perimeter of the geometry.
Definition: qgsabstractgeometry.cpp:137
QgsAbstractGeometry::isMeasure
bool isMeasure() const
Returns true if the geometry contains m values.
Definition: qgsabstractgeometry.h:211
QgsWkbTypes::isCurvedType
static bool isCurvedType(Type type)
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:880
QgsDistanceArea::computeSpheroidProject
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.
Definition: qgsdistancearea.cpp:396
QgsDistanceArea::latitudeGeodesicCrossesAntimeridian
double latitudeGeodesicCrossesAntimeridian(const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine) const
Calculates the latitude at which the geodesic line joining p1 and p2 crosses the antimeridian (longit...
Definition: qgsdistancearea.cpp:460
QgsEllipsoidUtils::EllipsoidParameters::inverseFlattening
double inverseFlattening
Inverse flattening.
Definition: qgsellipsoidutils.h:53
qgsgeometrycollection.h
qgsexception.h
QgsPointXY::setX
void setX(double x)
Sets the x value of the point.
Definition: qgspointxy.h:107
qgsdistancearea.h
QgsCurve::points
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
QgsDistanceArea::formatArea
static QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
Definition: qgsdistancearea.cpp:1137
qgslogger.h
QgsEllipsoidUtils::EllipsoidParameters::semiMajor
double semiMajor
Semi-major axis.
Definition: qgsellipsoidutils.h:45
QgsCurve::curveToLine
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.
qgscoordinatereferencesystem.h
QgsPoint::x
double x
Definition: qgspoint.h:58
qgsmultisurface.h
qgspointxy.h
qgsmultilinestring.h
qgsmessagelog.h
qgssurface.h
QgsCoordinateTransform::transformInPlace
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transforms an array of x, y and z double coordinates in place, from the source CRS to the destination...
Definition: qgscoordinatetransform.cpp:313
QgsDistanceArea::willUseEllipsoid
bool willUseEllipsoid() const
Returns whether calculations will use the ellipsoid.
Definition: qgsdistancearea.cpp:55