Expose method for calculating latitude geodesic crosses date line to …

…public QgsDistanceArea API
qgis · Jan 8, 2019 · adc5c17 · adc5c17
1 parent ebd1044
commit adc5c17
Show file tree

Hide file tree

Showing 3 changed files with 31 additions and 2 deletions.
diff --git a/python/core/auto_generated/qgsdistancearea.sip.in b/python/core/auto_generated/qgsdistancearea.sip.in
@@ -388,6 +388,19 @@ If the geodesic line crosses the international date line and ``breakLine`` is tr
 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
+
+    double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
+%Docstring
+Calculates the latitude at which the geodesic line joining ``p1`` and ``p2`` crosses
+the international date line (longitude +/- 180 degrees).
+
+The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
+
+``p1`` and ``p2`` must be in the ellipsoidCrs() of this QgsDistanceArea object. The returned latitude
+will also be in this same CRS.
+
 .. versionadded:: 3.6
 %End
 

diff --git a/src/core/qgsdistancearea.cpp b/src/core/qgsdistancearea.cpp
@@ -457,8 +457,10 @@ QgsPointXY QgsDistanceArea::computeSpheroidProject(
   return QgsPointXY( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
 }
 
-double calculateLatitudeGeodesicIntersectionAtDateLine( QgsPointXY p1, QgsPointXY p2, GeographicLib::Geodesic &geod )
+double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1, const QgsPointXY &pp2 ) const
 {
+  QgsPointXY p1 = pp1;
+  QgsPointXY p2 = pp2;
   if ( p1.x() < -120 )
     p1.setX( p1.x() + 360 );
   if ( p2.x() < -120 )
@@ -474,6 +476,7 @@ double calculateLatitudeGeodesicIntersectionAtDateLine( QgsPointXY p1, QgsPointX
   double lat = p2y;
   double lon = p2x;
 
+  GeographicLib::Geodesic geod( mSemiMajor, 1 / mInvFlattening );
   GeographicLib::GeodesicLine line = geod.InverseLine( p1y, p1x, p2y, p2x );
   double intersectionDist = line.Distance();
 
@@ -572,7 +575,7 @@ QList< QList<QgsPointXY> > QgsDistanceArea::geodesicLine( const QgsPointXY &p1,
       // 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.
-      double lat180 = calculateLatitudeGeodesicIntersectionAtDateLine( currentPart.constLast(), QgsPointXY( lon, lat ), geod );
+      double lat180 = latitudeGeodesicCrossesDateLine( currentPart.constLast(), QgsPointXY( lon, lat ) );
 
       try
       {

diff --git a/src/core/qgsdistancearea.h b/src/core/qgsdistancearea.h
@@ -319,6 +319,19 @@ class CORE_EXPORT QgsDistanceArea
      */
     QList< QList< QgsPointXY > > geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine = false ) const;
 
+    /**
+     * Calculates the latitude at which the geodesic line joining \a p1 and \a p2 crosses
+     * the international date line (longitude +/- 180 degrees).
+     *
+     * The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
+     *
+     * \a p1 and \a p2 must be in the ellipsoidCrs() of this QgsDistanceArea object. The returned latitude
+     * will also be in this same CRS.
+     *
+     * \since QGIS 3.6
+     */
+    double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
+
   private:
 
     /**