QGIS API Documentation  3.8.0-Zanzibar (11aff65)
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
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:486
QgsPolygon * surfaceToPolygon() const override
Gets a polygon representation of this surface.
Definition: qgspolygon.cpp:270
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:265
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.
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:111
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
Definition: qgspointxy.cpp:40
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:299
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:71
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 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:225
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:666
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
Definition: qgis.h:550
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:70
virtual double area() const
Returns the 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:190
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:104
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:213
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:224
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:54
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.
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:609
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 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:497
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:80
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