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