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