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