Skip to content

Commit

Permalink
Add method to calculate the geodesic line joining two points to QgsDi…
Browse files Browse the repository at this point in the history
…stanceArea

Using geographiclib to calculate the line
  • Loading branch information
nyalldawson committed Jan 8, 2019
1 parent 148505d commit fdea61a
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 0 deletions.
21 changes: 21 additions & 0 deletions python/core/auto_generated/qgsdistancearea.sip.in
Original file line number Diff line number Diff line change
Expand Up @@ -353,6 +353,27 @@ Datum of Australia Technical Manual", Chapter 4.
--> latitudes outside these bounds cause the calculations to become unstable and can return invalid results

.. versionadded:: 3.0
%End

QList< QList< QgsPointXY > > geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine = false ) const;
%Docstring
Calculates the geodesic line between ``p1`` and ``p2``, which represents the shortest path on the
ellipsoid between these two points.

The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.

``p1`` and ``p2`` must be in the sourceCrs() of this QgsDistanceArea object. The returned line
will also be in this same CRS.

The ``interval`` parameter gives the maximum distance between points on the computed line.
This argument is always specified in meters. A shorter distance results in a denser line,
at the cost of extra computing time.

If the geodesic line crosses the international date line and ``breakLine`` is true, then
the line will be split into two parts, broken at the date line. In this case the function
will return two lines, corresponding to the portions at either side of the date line.

.. versionadded:: 3.6
%End

};
Expand Down
63 changes: 63 additions & 0 deletions src/core/qgsdistancearea.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,9 @@
#include "qgsunittypes.h"
#include "qgsexception.h"

#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/GeodesicLine.hpp>

#define DEG2RAD(x) ((x)*M_PI/180)
#define RAD2DEG(r) (180.0 * (r) / M_PI)
#define POW2(x) ((x)*(x))
Expand Down Expand Up @@ -454,6 +457,66 @@ QgsPointXY QgsDistanceArea::computeSpheroidProject(
return QgsPointXY( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
}

QList< QList<QgsPointXY> > QgsDistanceArea::geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, const double interval, const bool breakLine ) const
{
if ( !willUseEllipsoid() )
{
return QList< QList< QgsPointXY > >() << ( QList< QgsPointXY >() << p1 << p2 );
}

GeographicLib::Geodesic geod( mSemiMajor, 1 / mInvFlattening );

QgsPointXY pp1, pp2;
try
{
pp1 = mCoordTransform.transform( p1 );
pp2 = mCoordTransform.transform( p2 );
}
catch ( QgsCsException & )
{
QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate geodesic line." ) );
return QList< QList< QgsPointXY > >();
}

GeographicLib::GeodesicLine line = geod.InverseLine( pp1.y(), pp1.x(), pp2.y(), pp2.x() );
const double totalDist = line.Distance();

QList< QList< QgsPointXY > > res;
QList< QgsPointXY > currentPart;
currentPart << p1;
double d = interval;
double prevLon = p1.x();
while ( d < totalDist )
{
double lat, lon;
( void )line.Position( d, lat, lon );

if ( breakLine && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
{
// TODO -- when breaking the geodesic at the date line, we need to calculate the latitude
// at which the geodesic intersects the date line, and add points to both line segments at this latitude
// on the date line. See Baselga & Martinez-Llario 2017 for an algorithm to calculate geodesic-geodesic intersections on ellipsoids.

res << currentPart;
currentPart.clear();
}

prevLon = lon;

try
{
currentPart << mCoordTransform.transform( QgsPointXY( lon, lat ), QgsCoordinateTransform::ReverseTransform );
}
catch ( QgsCsException & )
{
QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point." ) );
}
d += interval;
}
res << currentPart;
return res;
}

QgsUnitTypes::DistanceUnit QgsDistanceArea::lengthUnits() const
{
return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();
Expand Down
21 changes: 21 additions & 0 deletions src/core/qgsdistancearea.h
Original file line number Diff line number Diff line change
Expand Up @@ -288,6 +288,27 @@ class CORE_EXPORT QgsDistanceArea
*/
QgsPointXY computeSpheroidProject( const QgsPointXY &p1, double distance = 1, double azimuth = M_PI_2 ) const;

/**
* Calculates the geodesic line between \a p1 and \a p2, which represents the shortest path on the
* ellipsoid between these two points.
*
* The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
*
* \a p1 and \a p2 must be in the sourceCrs() of this QgsDistanceArea object. The returned line
* will also be in this same CRS.
*
* The \a interval parameter gives the maximum distance between points on the computed line.
* This argument is always specified in meters. A shorter distance results in a denser line,
* at the cost of extra computing time.
*
* If the geodesic line crosses the international date line and \a breakLine is true, then
* the line will be split into two parts, broken at the date line. In this case the function
* will return two lines, corresponding to the portions at either side of the date line.
*
* \since QGIS 3.6
*/
QList< QList< QgsPointXY > > geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine = false ) const;

private:

/**
Expand Down

0 comments on commit fdea61a

Please sign in to comment.