QGIS API Documentation  3.37.0-Master (a5b4d9743e8)
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 "qgslinestring.h"
31 #include "qgspolygon.h"
32 #include "qgssurface.h"
33 #include "qgsunittypes.h"
34 #include "qgsexception.h"
35 #include "qgsmultilinestring.h"
36 
37 #include <geodesic.h>
38 
39 #define DEG2RAD(x) ((x)*M_PI/180)
40 #define RAD2DEG(r) (180.0 * (r) / M_PI)
41 #define POW2(x) ((x)*(x))
42 
44 {
45  // init with default settings
46  mSemiMajor = -1.0;
47  mSemiMinor = -1.0;
48  mInvFlattening = -1.0;
49  const QgsCoordinateTransformContext context; // this is ok - by default we have a source/dest of WGS84, so no reprojection takes place
50  setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), context ); // WGS 84
51  setEllipsoid( geoNone() );
52 }
53 
55 
57  : mCoordTransform( other.mCoordTransform )
58  , mEllipsoid( other.mEllipsoid )
59  , mSemiMajor( other.mSemiMajor )
60  , mSemiMinor( other.mSemiMinor )
61  , mInvFlattening( other.mInvFlattening )
62 {
63  computeAreaInit();
64 }
65 
67 {
68  mCoordTransform = other.mCoordTransform;
69  mEllipsoid = other.mEllipsoid;
70  mSemiMajor = other.mSemiMajor;
71  mSemiMinor = other.mSemiMinor;
72  mInvFlattening = other.mInvFlattening;
73  computeAreaInit();
74  return *this;
75 }
76 
78 {
79  return mEllipsoid != geoNone();
80 }
81 
83 {
84  mCoordTransform.setContext( context );
85  mCoordTransform.setSourceCrs( srcCRS );
86 }
87 
88 bool QgsDistanceArea::setEllipsoid( const QString &ellipsoid )
89 {
90  // Shortcut if ellipsoid is none.
91  if ( ellipsoid == geoNone() )
92  {
93  mEllipsoid = geoNone();
94  mGeod.reset();
95  return true;
96  }
97 
99  if ( !params.valid )
100  {
101  mGeod.reset();
102  return false;
103  }
104  else
105  {
106  mEllipsoid = ellipsoid;
107  setFromParams( params );
108  return true;
109  }
110 }
111 
112 // Inverse flattening is calculated with invf = a/(a-b)
113 // Also, b = a-(a/invf)
114 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
115 {
116  mEllipsoid = QStringLiteral( "PARAMETER:%1:%2" ).arg( qgsDoubleToString( semiMajor ), qgsDoubleToString( semiMinor ) );
117  mSemiMajor = semiMajor;
118  mSemiMinor = semiMinor;
119  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
120 
121  computeAreaInit();
122 
123  return true;
124 }
125 
126 double QgsDistanceArea::measure( const QgsAbstractGeometry *geomV2, MeasureType type ) const
127 {
128  if ( !geomV2 )
129  {
130  return 0.0;
131  }
132 
133  const int geomDimension = geomV2->dimension();
134  if ( geomDimension <= 0 )
135  {
136  return 0.0;
137  }
138 
139  MeasureType measureType = type;
140  if ( measureType == Default )
141  {
142  measureType = ( geomDimension == 1 ? Length : Area );
143  }
144 
145  if ( !willUseEllipsoid() )
146  {
147  //no transform required
148  if ( measureType == Length )
149  {
150  return geomV2->length();
151  }
152  else
153  {
154  return geomV2->area();
155  }
156  }
157  else
158  {
159  //multigeom is sum of measured parts
160  const QgsGeometryCollection *collection = qgsgeometry_cast<const QgsGeometryCollection *>( geomV2 );
161  if ( collection )
162  {
163  double sum = 0;
164  for ( int i = 0; i < collection->numGeometries(); ++i )
165  {
166  sum += measure( collection->geometryN( i ), measureType );
167  }
168  return sum;
169  }
170 
171  if ( measureType == Length )
172  {
173  const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( geomV2 );
174  if ( !curve )
175  {
176  return 0.0;
177  }
178 
179  QgsLineString *lineString = curve->curveToLine();
180  const double length = measureLine( lineString );
181  delete lineString;
182  return length;
183  }
184  else
185  {
186  const QgsSurface *surface = qgsgeometry_cast<const QgsSurface *>( geomV2 );
187  if ( !surface )
188  return 0.0;
189 
190  QgsPolygon *polygon = surface->surfaceToPolygon();
191 
192  double area = 0;
193  const QgsCurve *outerRing = polygon->exteriorRing();
194  area += measurePolygon( outerRing );
195 
196  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
197  {
198  const QgsCurve *innerRing = polygon->interiorRing( i );
199  area -= measurePolygon( innerRing );
200  }
201  delete polygon;
202  return area;
203  }
204  }
205 }
206 
207 double QgsDistanceArea::measureArea( const QgsGeometry &geometry ) const
208 {
209  if ( geometry.isNull() )
210  return 0.0;
211 
212  const QgsAbstractGeometry *geomV2 = geometry.constGet();
213  return measure( geomV2, Area );
214 }
215 
216 double QgsDistanceArea::measureLength( const QgsGeometry &geometry ) const
217 {
218  if ( geometry.isNull() )
219  return 0.0;
220 
221  const QgsAbstractGeometry *geomV2 = geometry.constGet();
222  return measure( geomV2, Length );
223 }
224 
225 double QgsDistanceArea::measurePerimeter( const QgsGeometry &geometry ) const
226 {
227  if ( geometry.isNull() )
228  return 0.0;
229 
230  const QgsAbstractGeometry *geomV2 = geometry.constGet();
231  if ( !geomV2 || geomV2->dimension() < 2 )
232  {
233  return 0.0;
234  }
235 
236  if ( !willUseEllipsoid() )
237  {
238  return geomV2->perimeter();
239  }
240 
241  //create list with (single) surfaces
242  QVector< const QgsSurface * > surfaces;
243  const QgsSurface *surf = qgsgeometry_cast<const QgsSurface *>( geomV2 );
244  if ( surf )
245  {
246  surfaces.append( surf );
247  }
248  const QgsMultiSurface *multiSurf = qgsgeometry_cast<const QgsMultiSurface *>( geomV2 );
249  if ( multiSurf )
250  {
251  surfaces.reserve( ( surf ? 1 : 0 ) + multiSurf->numGeometries() );
252  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
253  {
254  surfaces.append( static_cast<const QgsSurface *>( multiSurf->geometryN( i ) ) );
255  }
256  }
257 
258  double length = 0;
259  QVector<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
260  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
261  {
262  if ( !*surfaceIt )
263  {
264  continue;
265  }
266 
267  QgsPolygon *poly = ( *surfaceIt )->surfaceToPolygon();
268  const QgsCurve *outerRing = poly->exteriorRing();
269  if ( outerRing )
270  {
271  length += measure( outerRing );
272  }
273  const int nInnerRings = poly->numInteriorRings();
274  for ( int i = 0; i < nInnerRings; ++i )
275  {
276  length += measure( poly->interiorRing( i ) );
277  }
278  delete poly;
279  }
280  return length;
281 }
282 
283 double QgsDistanceArea::measureLine( const QgsCurve *curve ) const
284 {
285  if ( !curve )
286  {
287  return 0.0;
288  }
289 
290  QgsPointSequence linePointsV2;
291  QVector<QgsPointXY> linePoints;
292  curve->points( linePointsV2 );
293  QgsGeometry::convertPointList( linePointsV2, linePoints );
294  return measureLine( linePoints );
295 }
296 
297 double QgsDistanceArea::measureLine( const QVector<QgsPointXY> &points ) const
298 {
299  if ( points.size() < 2 )
300  return 0;
301 
302  double total = 0;
303  QgsPointXY p1, p2;
304 
305  if ( willUseEllipsoid() )
306  {
307  if ( !mGeod )
308  computeAreaInit();
309  Q_ASSERT_X( static_cast<bool>( mGeod ), "QgsDistanceArea::measureLine()", "Error creating geod_geodesic object" );
310  if ( !mGeod )
311  return 0;
312  }
313 
314  try
315  {
316  if ( willUseEllipsoid() )
317  p1 = mCoordTransform.transform( points[0] );
318  else
319  p1 = points[0];
320 
321  for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
322  {
323  if ( willUseEllipsoid() )
324  {
325  p2 = mCoordTransform.transform( *i );
326 
327  double distance = 0;
328  double azimuth1 = 0;
329  double azimuth2 = 0;
330  geod_inverse( mGeod.get(), p1.y(), p1.x(), p2.y(), p2.x(), &distance, &azimuth1, &azimuth2 );
331  total += distance;
332  }
333  else
334  {
335  p2 = *i;
336  total += measureLine( p1, p2 );
337  }
338 
339  p1 = p2;
340  }
341 
342  return total;
343  }
344  catch ( QgsCsException &cse )
345  {
346  Q_UNUSED( cse )
347  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
348  return 0.0;
349  }
350 
351 }
352 
353 double QgsDistanceArea::measureLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const
354 {
355  double result;
356 
357  if ( willUseEllipsoid() )
358  {
359  if ( !mGeod )
360  computeAreaInit();
361  Q_ASSERT_X( static_cast<bool>( mGeod ), "QgsDistanceArea::measureLine()", "Error creating geod_geodesic object" );
362  if ( !mGeod )
363  return 0;
364  }
365 
366  try
367  {
368  QgsPointXY pp1 = p1, pp2 = p2;
369 
370  QgsDebugMsgLevel( QStringLiteral( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
371  if ( willUseEllipsoid() )
372  {
373  QgsDebugMsgLevel( QStringLiteral( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
374  QgsDebugMsgLevel( QStringLiteral( "From proj4 : %1" ).arg( mCoordTransform.sourceCrs().toProj() ), 4 );
375  QgsDebugMsgLevel( QStringLiteral( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj() ), 4 );
376  pp1 = mCoordTransform.transform( p1 );
377  pp2 = mCoordTransform.transform( p2 );
378  QgsDebugMsgLevel( QStringLiteral( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
379 
380  double azimuth1 = 0;
381  double azimuth2 = 0;
382  geod_inverse( mGeod.get(), pp1.y(), pp1.x(), pp2.y(), pp2.x(), &result, &azimuth1, &azimuth2 );
383  }
384  else
385  {
386  QgsDebugMsgLevel( QStringLiteral( "Cartesian calculation on canvas coordinates" ), 4 );
387  result = p2.distance( p1 );
388  }
389  }
390  catch ( QgsCsException &cse )
391  {
392  Q_UNUSED( cse )
393  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
394  result = 0.0;
395  }
396  QgsDebugMsgLevel( QStringLiteral( "The result was %1" ).arg( result ), 3 );
397  return result;
398 }
399 
400 double QgsDistanceArea::measureLineProjected( const QgsPointXY &p1, double distance, double azimuth, QgsPointXY *projectedPoint ) const
401 {
402  double result = 0.0;
403  QgsPointXY p2;
404  if ( mCoordTransform.sourceCrs().isGeographic() && willUseEllipsoid() )
405  {
406  p2 = computeSpheroidProject( p1, distance, azimuth );
407  result = p1.distance( p2 );
408  }
409  else // Cartesian coordinates
410  {
411  result = distance; // Avoid rounding errors when using meters [return as sent]
412  if ( sourceCrs().mapUnits() != Qgis::DistanceUnit::Meters )
413  {
414  distance = ( distance * QgsUnitTypes::fromUnitToUnitFactor( Qgis::DistanceUnit::Meters, sourceCrs().mapUnits() ) );
415  result = p1.distance( p2 );
416  }
417  p2 = p1.project( distance, azimuth );
418  }
419  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]" )
420  .arg( QString::number( distance, 'f', 7 ),
422  QString::number( result, 'f', 7 ),
423  mCoordTransform.sourceCrs().isGeographic() ? QStringLiteral( "Geographic" ) : QStringLiteral( "Cartesian" ),
424  QgsUnitTypes::toString( sourceCrs().mapUnits() ) )
425  .arg( azimuth )
426  .arg( p1.asWkt(),
427  p2.asWkt(),
428  sourceCrs().description(),
429  mEllipsoid )
430  .arg( sourceCrs().isGeographic() )
431  .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 );
432  if ( projectedPoint )
433  {
434  *projectedPoint = QgsPointXY( p2 );
435  }
436  return result;
437 }
438 
440  const QgsPointXY &p1, double distance, double azimuth ) const
441 {
442  if ( !mGeod )
443  computeAreaInit();
444  if ( !mGeod )
445  return QgsPointXY();
446 
447  double lat2 = 0;
448  double lon2 = 0;
449  double azimuth2 = 0;
450  geod_direct( mGeod.get(), p1.y(), p1.x(), RAD2DEG( azimuth ), distance, &lat2, &lon2, &azimuth2 );
451 
452  return QgsPointXY( lon2, lat2 );
453 }
454 
455 double QgsDistanceArea::latitudeGeodesicCrossesAntimeridian( const QgsPointXY &pp1, const QgsPointXY &pp2, double &fractionAlongLine ) const
456 {
457  QgsPointXY p1 = pp1;
458  QgsPointXY p2 = pp2;
459  if ( p1.x() < -120 )
460  p1.setX( p1.x() + 360 );
461  if ( p2.x() < -120 )
462  p2.setX( p2.x() + 360 );
463 
464  // we need p2.x() > 180 and p1.x() < 180
465  double p1x = p1.x() < 180 ? p1.x() : p2.x();
466  double p1y = p1.x() < 180 ? p1.y() : p2.y();
467  double p2x = p1.x() < 180 ? p2.x() : p1.x();
468  double p2y = p1.x() < 180 ? p2.y() : p1.y();
469  // lat/lon are our candidate intersection position - we want this to get as close to 180 as possible
470  // the first candidate is p2
471  double lat = p2y;
472  double lon = p2x;
473 
474  if ( mEllipsoid == geoNone() )
475  {
476  fractionAlongLine = ( 180 - p1x ) / ( p2x - p1x );
477  if ( p1.x() >= 180 )
478  fractionAlongLine = 1 - fractionAlongLine;
479  return p1y + ( 180 - p1x ) / ( p2x - p1x ) * ( p2y - p1y );
480  }
481 
482  if ( !mGeod )
483  computeAreaInit();
484  Q_ASSERT_X( static_cast<bool>( mGeod ), "QgsDistanceArea::latitudeGeodesicCrossesAntimeridian()", "Error creating geod_geodesic object" );
485  if ( !mGeod )
486  return 0;
487 
488  geod_geodesicline line;
489  geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
490 
491  const double totalDist = line.s13;
492  double intersectionDist = line.s13;
493 
494  int iterations = 0;
495  double t = 0;
496  // iterate until our intersection candidate is within ~1 mm of the antimeridian (or too many iterations happened)
497  while ( std::fabs( lon - 180.0 ) > 0.00000001 && iterations < 100 )
498  {
499  if ( iterations > 0 && std::fabs( p2x - p1x ) > 5 )
500  {
501  // if we have too large a range of longitudes, we use a binary search to narrow the window -- this ensures we will converge
502  if ( lon < 180 )
503  {
504  p1x = lon;
505  p1y = lat;
506  }
507  else
508  {
509  p2x = lon;
510  p2y = lat;
511  }
512  QgsDebugMsgLevel( QStringLiteral( "Narrowed window to %1, %2 - %3, %4" ).arg( p1x ).arg( p1y ).arg( p2x ).arg( p2y ), 4 );
513 
514  geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
515  intersectionDist = line.s13 * 0.5;
516  }
517  else
518  {
519  // we have a sufficiently narrow window -- use Newton's method
520  // adjust intersection distance by fraction of how close the previous candidate was to 180 degrees longitude -
521  // this helps us close in to the correct longitude quickly
522  intersectionDist *= ( 180.0 - p1x ) / ( lon - p1x );
523  }
524 
525  // now work out the point on the geodesic this far from p1 - this becomes our new candidate for crossing the antimeridian
526 
527  geod_position( &line, intersectionDist, &lat, &lon, &t );
528  // we don't want to wrap longitudes > 180 around)
529  if ( lon < 0 )
530  lon += 360;
531 
532  iterations++;
533  QgsDebugMsgLevel( QStringLiteral( "After %1 iterations lon is %2, lat is %3, dist from p1: %4" ).arg( iterations ).arg( lon ).arg( lat ).arg( intersectionDist ), 4 );
534  }
535 
536  fractionAlongLine = intersectionDist / totalDist;
537  if ( p1.x() >= 180 )
538  fractionAlongLine = 1 - fractionAlongLine;
539 
540  // either converged on 180 longitude or hit too many iterations
541  return lat;
542 }
543 
545 {
547  return geometry;
548 
549  QgsGeometry g = geometry;
550  // TODO - avoid segmentization of curved geometries (if this is even possible!)
551  if ( QgsWkbTypes::isCurvedType( g.wkbType() ) )
553 
554  std::unique_ptr< QgsMultiLineString > res = std::make_unique< QgsMultiLineString >();
555  for ( auto part = g.const_parts_begin(); part != g.const_parts_end(); ++part )
556  {
557  const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( *part );
558  if ( !line )
559  continue;
560  if ( line->isEmpty() )
561  {
562  continue;
563  }
564 
565  const std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >();
566  try
567  {
568  double x = 0;
569  double y = 0;
570  double z = 0;
571  double m = 0;
572  QVector< QgsPoint > newPoints;
573  newPoints.reserve( line->numPoints() );
574  double prevLon = 0;
575  double prevLat = 0;
576  double lon = 0;
577  double lat = 0;
578  double prevZ = 0;
579  double prevM = 0;
580  for ( int i = 0; i < line->numPoints(); i++ )
581  {
582  QgsPoint p = line->pointN( i );
583  x = p.x();
584  if ( mCoordTransform.sourceCrs().isGeographic() )
585  {
586  x = std::fmod( x, 360.0 );
587  if ( x > 180 )
588  x -= 360;
589  p.setX( x );
590  }
591  y = p.y();
592  lon = x;
593  lat = y;
594  mCoordTransform.transformInPlace( lon, lat, z );
595 
596  //test if we crossed the antimeridian in this segment
597  if ( i > 0 && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
598  {
599  // we did!
600  // when crossing the antimeridian, we need to calculate the latitude
601  // at which the geodesic intersects the antimeridian
602  double fract = 0;
603  const double lat180 = latitudeGeodesicCrossesAntimeridian( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ), fract );
604  if ( line->is3D() )
605  {
606  z = prevZ + ( p.z() - prevZ ) * fract;
607  }
608  if ( line->isMeasure() )
609  {
610  m = prevM + ( p.m() - prevM ) * fract;
611  }
612 
613  QgsPointXY antiMeridianPoint;
614  if ( prevLon < -120 )
615  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), Qgis::TransformDirection::Reverse );
616  else
617  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), Qgis::TransformDirection::Reverse );
618 
619  QgsPoint newPoint( antiMeridianPoint );
620  if ( line->is3D() )
621  newPoint.addZValue( z );
622  if ( line->isMeasure() )
623  newPoint.addMValue( m );
624 
625  if ( std::isfinite( newPoint.x() ) && std::isfinite( newPoint.y() ) )
626  {
627  newPoints << newPoint;
628  }
629  res->addGeometry( new QgsLineString( newPoints ) );
630 
631  newPoints.clear();
632  newPoints.reserve( line->numPoints() - i + 1 );
633 
634  if ( lon < -120 )
635  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( -180, lat180 ), Qgis::TransformDirection::Reverse );
636  else
637  antiMeridianPoint = mCoordTransform.transform( QgsPointXY( 180, lat180 ), Qgis::TransformDirection::Reverse );
638 
639  if ( std::isfinite( antiMeridianPoint.x() ) && std::isfinite( antiMeridianPoint.y() ) )
640  {
641  // we want to keep the previously calculated z/m value for newPoint, if present. They're the same each
642  // of the antimeridian split
643  newPoint.setX( antiMeridianPoint.x() );
644  newPoint.setY( antiMeridianPoint.y() );
645  newPoints << newPoint;
646  }
647  }
648  newPoints << p;
649 
650  prevLon = lon;
651  prevLat = lat;
652  if ( line->is3D() )
653  prevZ = p.z();
654  if ( line->isMeasure() )
655  prevM = p.m();
656  }
657  res->addGeometry( new QgsLineString( newPoints ) );
658  }
659  catch ( QgsCsException & )
660  {
661  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform linestring. Unable to calculate break point." ) );
662  res->addGeometry( line->clone() );
663  break;
664  }
665  }
666 
667  return QgsGeometry( std::move( res ) );
668 }
669 
670 
671 QVector< QVector<QgsPointXY> > QgsDistanceArea::geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, const double interval, const bool breakLine ) const
672 {
673  if ( !willUseEllipsoid() )
674  {
675  return QVector< QVector< QgsPointXY > >() << ( QVector< QgsPointXY >() << p1 << p2 );
676  }
677 
678  if ( !mGeod )
679  computeAreaInit();
680  if ( !mGeod )
681  return QVector< QVector< QgsPointXY > >();
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, mGeod.get(), 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  const 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 ), Qgis::TransformDirection::Reverse );
735  else
736  p = mCoordTransform.transform( QgsPointXY( 180, lat180 ), Qgis::TransformDirection::Reverse );
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 ), Qgis::TransformDirection::Reverse );
753  else
754  p = mCoordTransform.transform( QgsPointXY( 180, lat180 ), Qgis::TransformDirection::Reverse );
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 ), Qgis::TransformDirection::Reverse );
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() ? Qgis::DistanceUnit::Meters : 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 
852  if ( !mGeod )
853  computeAreaInit();
854  Q_ASSERT_X( static_cast<bool>( mGeod ), "QgsDistanceArea::bearing()", "Error creating geod_geodesic object" );
855  if ( !mGeod )
856  return 0;
857 
858  double distance = 0;
859  double azimuth1 = 0;
860  double azimuth2 = 0;
861  geod_inverse( mGeod.get(), pp1.y(), pp1.x(), pp2.y(), pp2.x(), &distance, &azimuth1, &azimuth2 );
862 
863  bearing = DEG2RAD( azimuth1 );
864  }
865  else //compute simple planar azimuth
866  {
867  const double dx = p2.x() - p1.x();
868  const double dy = p2.y() - p1.y();
869  // Note: the prototype of std::atan2 is (y,x), to return the angle of
870  // vector (x,y) from the horizontal axis in counter-clock-wise orientation.
871  // But a bearing is expressed in clock-wise order from the vertical axis, so
872  // M_PI / 2 - std::atan2( dy, dx ) == std::atan2( dx, dy )
873  bearing = std::atan2( dx, dy );
874  }
875 
876  return bearing;
877 }
878 
879 void QgsDistanceArea::computeAreaInit() const
880 {
881  //don't try to perform calculations if no ellipsoid
882  if ( mEllipsoid == geoNone() )
883  {
884  mGeod.reset();
885  return;
886  }
887 
888  mGeod.reset( new geod_geodesic() );
889  geod_init( mGeod.get(), mSemiMajor, 1 / mInvFlattening );
890 }
891 
892 void QgsDistanceArea::setFromParams( const QgsEllipsoidUtils::EllipsoidParameters &params )
893 {
894  if ( params.useCustomParameters )
895  {
896  setEllipsoid( params.semiMajor, params.semiMinor );
897  }
898  else
899  {
900  mSemiMajor = params.semiMajor;
901  mSemiMinor = params.semiMinor;
902  mInvFlattening = params.inverseFlattening;
903  mCoordTransform.setDestinationCrs( params.crs );
904  computeAreaInit();
905  }
906 }
907 
908 double QgsDistanceArea::computePolygonArea( const QVector<QgsPointXY> &points ) const
909 {
910  if ( points.isEmpty() )
911  {
912  return 0;
913  }
914 
915  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
916  if ( !willUseEllipsoid() )
917  {
918  return computePolygonFlatArea( points );
919  }
920 
921  if ( !mGeod )
922  computeAreaInit();
923  Q_ASSERT_X( static_cast<bool>( mGeod ), "QgsDistanceArea::computePolygonArea()", "Error creating geod_geodesic object" );
924  if ( !mGeod )
925  return 0;
926 
927  struct geod_polygon p;
928  geod_polygon_init( &p, 0 );
929 
930  const bool isClosed = points.constFirst() == points.constLast();
931 
932  /* GeographicLib does not need a closed ring,
933  * see example for geod_polygonarea() in geodesic.h */
934  /* add points in reverse order */
935  int i = points.size();
936  while ( ( isClosed && --i ) || ( !isClosed && --i >= 0 ) )
937  geod_polygon_addpoint( mGeod.get(), &p, points.at( i ).y(), points.at( i ).x() );
938 
939  double area = 0;
940  double perimeter = 0;
941  geod_polygon_compute( mGeod.get(), &p, 0, 1, &area, &perimeter );
942 
943  return std::fabs( area );
944 }
945 
946 double QgsDistanceArea::computePolygonFlatArea( const QVector<QgsPointXY> &points ) const
947 {
948  // Normal plane area calculations.
949  double area = 0.0;
950  int i, size;
951 
952  size = points.size();
953 
954  // QgsDebugMsgLevel("New area calc, nr of points: " + QString::number(size), 2);
955  for ( i = 0; i < size; i++ )
956  {
957  // QgsDebugMsgLevel("Area from point: " + (points[i]).toString(2), 2);
958  // Using '% size', so that we always end with the starting point
959  // and thus close the polygon.
960  area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
961  }
962  // QgsDebugMsgLevel("Area from point: " + (points[i % size]).toString(2), 2);
963  area = area / 2.0;
964  return std::fabs( area ); // All areas are positive!
965 }
966 
967 QString QgsDistanceArea::formatDistance( double distance, int decimals, Qgis::DistanceUnit unit, bool keepBaseUnit )
968 {
969  return QgsUnitTypes::formatDistance( distance, decimals, unit, keepBaseUnit );
970 }
971 
972 QString QgsDistanceArea::formatArea( double area, int decimals, Qgis::AreaUnit unit, bool keepBaseUnit )
973 {
974  return QgsUnitTypes::formatArea( area, decimals, unit, keepBaseUnit );
975 }
976 
977 double QgsDistanceArea::convertLengthMeasurement( double length, Qgis::DistanceUnit toUnits ) const
978 {
979  // get the conversion factor between the specified units
980  const Qgis::DistanceUnit measureUnits = lengthUnits();
981  const double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
982 
983  const double result = length * factorUnits;
984  QgsDebugMsgLevel( QStringLiteral( "Converted length of %1 %2 to %3 %4" ).arg( length )
985  .arg( QgsUnitTypes::toString( measureUnits ) )
986  .arg( result )
987  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
988  return result;
989 }
990 
991 double QgsDistanceArea::convertAreaMeasurement( double area, Qgis::AreaUnit toUnits ) const
992 {
993  // get the conversion factor between the specified units
994  const Qgis::AreaUnit measureUnits = areaUnits();
995  const double factorUnits = QgsUnitTypes::fromUnitToUnitFactor( measureUnits, toUnits );
996 
997  const double result = area * factorUnits;
998  QgsDebugMsgLevel( QStringLiteral( "Converted area of %1 %2 to %3 %4" ).arg( area )
999  .arg( QgsUnitTypes::toString( measureUnits ) )
1000  .arg( result )
1001  .arg( QgsUnitTypes::toString( toUnits ) ), 3 );
1002  return result;
1003 }
DistanceUnit
Units of distance.
Definition: qgis.h:4090
AreaUnit
Units of area.
Definition: qgis.h:4128
@ SquareMeters
Square meters.
@ Reverse
Reverse/inverse transform (from destination to source)
Abstract base class for all geometries.
bool isMeasure() const
Returns true if the geometry contains m values.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual double perimeter() const
Returns the planar, 2-dimensional perimeter of the geometry.
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual double area() const
Returns the planar, 2-dimensional area of the geometry.
This class represents a coordinate reference system (CRS).
Q_GADGET Qgis::DistanceUnit mapUnits
QString toProj() const
Returns a Proj string representation of this CRS.
Contains information about the context in which a coordinate transform is executed.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from.
void setContext(const QgsCoordinateTransformContext &context)
Sets the context in which the coordinate transform should be calculated.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
QgsPointXY transform(const QgsPointXY &point, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const
Transform the point from the source CRS to the destination CRS.
void transformInPlace(double &x, double &y, double &z, Qgis::TransformDirection direction=Qgis::TransformDirection::Forward) const
Transforms an array of x, y and z double coordinates in place, from the source CRS to the destination...
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:67
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
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.
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
QgsDistanceArea & operator=(const QgsDistanceArea &other)
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...
static QString formatDistance(double distance, int decimals, Qgis::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
double convertLengthMeasurement(double length, Qgis::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
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...
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
QString ellipsoid() const
Returns ellipsoid's acronym.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
QgsGeometry splitGeometryAtAntimeridian(const QgsGeometry &geometry) const
Splits a (Multi)LineString geometry at the antimeridian (longitude +/- 180 degrees).
Qgis::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
double measurePolygon(const QVector< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
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...
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
QgsDistanceArea()
Constructor.
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.
double convertAreaMeasurement(double area, Qgis::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
static QString formatArea(double area, int decimals, Qgis::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
Qgis::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
Geometry collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
int numGeometries() const
Returns the number of geometries within the collection.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:162
QgsAbstractGeometry::const_part_iterator const_parts_begin() const
Returns STL-style const iterator pointing to the first part of the geometry.
Q_GADGET bool isNull
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
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.
QgsAbstractGeometry::const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary part after the last part of the geometry.
static void convertPointList(const QVector< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPoint.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
bool isEmpty() const override
Returns true if the geometry is empty.
int numPoints() const override
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Multi surface geometry collection.
A class to represent a 2D point.
Definition: qgspointxy.h:60
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
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
Definition: qgspointxy.cpp:51
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:207
QString asWkt() const
Returns the well known text representation for the point (e.g.
Definition: qgspointxy.cpp:69
double y
Definition: qgspointxy.h:64
Q_GADGET double x
Definition: qgspointxy.h:63
void setX(double x)
Sets the x value of the point.
Definition: qgspointxy.h:120
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
void setY(double y)
Sets the point's y-coordinate.
Definition: qgspoint.h:343
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:564
void clear() override
Clears the geometry, ie reset it to a null geometry.
Definition: qgspoint.cpp:356
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:553
void setX(double x)
Sets the point's x-coordinate.
Definition: qgspoint.h:332
Q_GADGET double x
Definition: qgspoint.h:52
double z
Definition: qgspoint.h:54
double m
Definition: qgspoint.h:55
double y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:33
QgsPolygon * surfaceToPolygon() const override
Gets a polygon representation of this surface.
Definition: qgspolygon.cpp:355
Surface geometry type.
Definition: qgssurface.h:34
virtual QgsPolygon * surfaceToPolygon() const =0
Gets a polygon representation of this surface.
static Q_INVOKABLE QString toString(Qgis::DistanceUnit unit)
Returns a translated string representing a distance unit.
static Q_INVOKABLE QString formatArea(double area, int decimals, Qgis::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
static Q_INVOKABLE double fromUnitToUnitFactor(Qgis::DistanceUnit fromUnit, Qgis::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
static Q_INVOKABLE QString formatDistance(double distance, int decimals, Qgis::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
static Q_INVOKABLE Qgis::AreaUnit distanceToAreaUnit(Qgis::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:862
static bool isCurvedType(Qgis::WkbType type)
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:806
CONSTLATIN1STRING geoNone()
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.h:5659
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:5089
QVector< QgsPoint > QgsPointSequence
#define RAD2DEG(r)
#define DEG2RAD(x)
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
Contains parameters for an ellipsoid.
bool valid
Whether ellipsoid parameters are valid.
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
double inverseFlattening
Inverse flattening.
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)