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