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