QGIS API Documentation  3.10.0-A Coruña (6c816b4204)
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
53 }
54 
56 {
57  return mEllipsoid != GEO_NONE;
58 }
59 
61 {
62  mCoordTransform.setContext( context );
63  mCoordTransform.setSourceCrs( srcCRS );
64 }
65 
67 {
68  // Shortcut if ellipsoid is none.
69  if ( ellipsoid == GEO_NONE )
70  {
71  mEllipsoid = GEO_NONE;
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().toProj4() ), 4 );
328  QgsDebugMsgLevel( QStringLiteral( "To proj4 : %1" ).arg( mCoordTransform.destinationCrs().toProj4() ), 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(),
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 == GEO_NONE )
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 == GEO_NONE )
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 }
bool isMeasure() const
Returns true if the geometry contains m values.
QgsGeometry splitGeometryAtAntimeridian(const QgsGeometry &geometry) const
Splits a (Multi)LineString geometry at the antimeridian (longitude +/- 180 degrees).
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
bool isEmpty() const override
Returns true if the geometry is empty.
double y
Definition: qgspoint.h:42
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)
static Q_INVOKABLE QString toString(QgsUnitTypes::DistanceUnit unit)
Returns a translated string representing a distance unit.
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:501
QgsPolygon * surfaceToPolygon() const override
Gets a polygon representation of this surface.
Definition: qgspolygon.cpp:271
void setContext(const QgsCoordinateTransformContext &context)
Sets the context in which the coordinate transform should be calculated.
QgsPointXY transform(const QgsPointXY &point, TransformDirection direction=ForwardTransform) const SIP_THROW(QgsCsException)
Transform the point from the source CRS to the destination CRS.
QgsWkbTypes::Type wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
double y
Definition: qgspointxy.h:48
A class to represent a 2D point.
Definition: qgspointxy.h:43
QgsAbstractGeometry::const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary part after the last part of the geometry...
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition: qgis.h:280
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
QgsAbstractGeometry::const_part_iterator const_parts_begin() const
Returns STL-style const iterator pointing to the first part of the geometry.
#define DEG2RAD(x)
QString toProj4() const
Returns a Proj4 string representation of this CRS.
QgsCoordinateReferenceSystem destinationCrs() const
Returns the destination coordinate reference system, which the transform will transform coordinates t...
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:122
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
Definition: qgspointxy.cpp:51
double convertLengthMeasurement(double length, QgsUnitTypes::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
double convertAreaMeasurement(double area, QgsUnitTypes::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
Contains parameters for an ellipsoid.
double measurePolygon(const QVector< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
void clear() override
Clears the geometry, ie reset it to a null geometry.
Definition: qgspoint.cpp:314
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. ...
Multi surface geometry collection.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:72
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...
int numPoints() const override
Returns the number of points in the curve.
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
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.
virtual QgsPolygon * surfaceToPolygon() const =0
Gets a polygon representation of this surface.
static Q_INVOKABLE QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source coordinate reference system, which the transform will transform coordinates from...
void setDestinationCrs(const QgsCoordinateReferenceSystem &crs)
Sets the destination coordinate reference system.
bool valid
Whether ellipsoid parameters are valid.
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:39
Geometry collection.
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
Definition: qgis.h:240
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).
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:812
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
Definition: qgis.h:575
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
QString description() const
Returns the descriptive name of the CRS, e.g., "WGS 84" or "GDA 94 / Vicgrid94".
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
virtual double area() const
Returns the planar, 2-dimensional area of the geometry.
Abstract base class for curved geometry type.
Definition: qgscurve.h:35
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
Abstract base class for all geometries.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
Contains information about the context in which a coordinate transform is executed.
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Definition: qgspointxy.h:196
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:37
#define POW2(x)
void setX(double x)
Sets the x value of the point.
Definition: qgspointxy.h:107
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
double x
Definition: qgspointxy.h:47
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
int numGeometries() const
Returns the number of geometries within the collection.
void setX(double x)
Sets the point&#39;s x-coordinate.
Definition: qgspoint.h:262
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
void setY(double y)
Sets the point&#39;s y-coordinate.
Definition: qgspoint.h:273
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...
DistanceUnit
Units of distance.
Definition: qgsunittypes.h:66
QVector< QgsPoint > QgsPointSequence
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
static void convertPointList(const QVector< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPoint.
QgsUnitTypes::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
void reserve(int size)
Attempts to allocate memory for at least size geometries.
QString ellipsoid() const
Returns ellipsoid&#39;s acronym.
Transform from destination to source CRS.
static Q_INVOKABLE QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
static bool isCurvedType(Type type)
Returns true if the WKB type is a curved type or can contain curved geometries.
Definition: qgswkbtypes.h:755
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...
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:43
virtual double perimeter() const
Returns the planar, 2-dimensional perimeter of the geometry.
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
This class represents a coordinate reference system (CRS).
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
double inverseFlattening
Inverse flattening.
static QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
Definition: qgspoint.cpp:512
double z
Definition: qgspoint.h:43
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
Custom exception class for Coordinate Reference System related exceptions.
Definition: qgsexception.h:65
Polygon geometry type.
Definition: qgspolygon.h:31
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...
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
const QgsCurve * exteriorRing() const
Returns the curve polygon&#39;s exterior ring.
static Q_INVOKABLE double fromUnitToUnitFactor(QgsUnitTypes::DistanceUnit fromUnit, QgsUnitTypes::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
QgsDistanceArea()
Constructor.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs)
Sets the source coordinate reference system.
static Q_INVOKABLE QgsUnitTypes::AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
AreaUnit
Units of area.
Definition: qgsunittypes.h:92
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
static QgsCoordinateReferenceSystem fromSrsId(long srsId)
Creates a CRS from a specified QGIS SRS ID.
double m
Definition: qgspoint.h:44
#define RAD2DEG(r)
QString asWkt() const
Returns the well known text representation for the point (e.g.
Definition: qgspointxy.cpp:69
static QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
double x
Definition: qgspoint.h:41