Add method to calculate the geodesic line joining two points to QgsDistanceArea

nyalldawson · nyalldawson · commit fdea61a97ced · 2019-01-08T12:32:50.000+10:00
Using geographiclib to calculate the line
diff --git a/python/core/auto_generated/qgsdistancearea.sip.in b/python/core/auto_generated/qgsdistancearea.sip.in
@@ -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
 
 };
diff --git a/src/core/qgsdistancearea.cpp b/src/core/qgsdistancearea.cpp
@@ -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))
@@ -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();
diff --git a/src/core/qgsdistancearea.h b/src/core/qgsdistancearea.h
@@ -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:
 
     /**
