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