Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
latitudeGeodesicCrossesDateLine also returns fraction of segment
along which the crossing occurs
  • Loading branch information
nyalldawson committed Jan 8, 2019
1 parent e286c09 commit 13ac907
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 41 deletions.
8 changes: 7 additions & 1 deletion python/core/auto_generated/qgsdistancearea.sip.in
Expand Up @@ -391,7 +391,7 @@ will return two lines, corresponding to the portions at either side of the date
.. versionadded:: 3.6
%End

double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine /Out/ ) const;
%Docstring
Calculates the latitude at which the geodesic line joining ``p1`` and ``p2`` crosses
the international date line (longitude +/- 180 degrees).
Expand All @@ -401,6 +401,12 @@ The ellipsoid settings defined on this QgsDistanceArea object will be used durin
``p1`` and ``p2`` must be in the ellipsoidCrs() of this QgsDistanceArea object. The returned latitude
will also be in this same CRS.

:param p1: Starting point, in ellipsoidCrs()
:param p2: Ending point, in ellipsoidCrs()

:return: - the latitude at which the geodesic crosses the date line
- fractionAlongLine: will be set to the fraction along the geodesic line joining ``p1`` to ``p2`` at which the date line crossing occurs.

.. versionadded:: 3.6
%End

Expand Down
8 changes: 6 additions & 2 deletions src/core/qgsdistancearea.cpp
Expand Up @@ -456,7 +456,7 @@ QgsPointXY QgsDistanceArea::computeSpheroidProject(
return QgsPointXY( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
}

double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1, const QgsPointXY &pp2 ) const
double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1, const QgsPointXY &pp2, double &fractionAlongLine ) const
{
QgsPointXY p1 = pp1;
QgsPointXY p2 = pp2;
Expand All @@ -481,6 +481,7 @@ double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1,
geod_geodesicline line;
geod_inverseline( &line, &geod, p1y, p1x, p2y, p2x, GEOD_ALL );

const double totalDist = line.s13;
double intersectionDist = line.s13;

int iterations = 0;
Expand Down Expand Up @@ -525,6 +526,8 @@ double QgsDistanceArea::latitudeGeodesicCrossesDateLine( const QgsPointXY &pp1,
QgsDebugMsgLevel( QStringLiteral( "After %1 iterations lon is %2, lat is %3, dist from p1: %4" ).arg( iterations ).arg( lon ).arg( lat ).arg( intersectionDist ), 4 );
}

fractionAlongLine = intersectionDist / totalDist;

// either converged on 180 longitude or hit too many iterations
return lat;
}
Expand Down Expand Up @@ -583,7 +586,8 @@ 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 = latitudeGeodesicCrossesDateLine( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ) );
double fraction;
double lat180 = latitudeGeodesicCrossesDateLine( QgsPointXY( prevLon, prevLat ), QgsPointXY( lon, lat ), fraction );

try
{
Expand Down
7 changes: 6 additions & 1 deletion src/core/qgsdistancearea.h
Expand Up @@ -328,9 +328,14 @@ class CORE_EXPORT QgsDistanceArea
* \a p1 and \a p2 must be in the ellipsoidCrs() of this QgsDistanceArea object. The returned latitude
* will also be in this same CRS.
*
* \param p1 Starting point, in ellipsoidCrs()
* \param p2 Ending point, in ellipsoidCrs()
* \param fractionAlongLine will be set to the fraction along the geodesic line joining \a p1 to \a p2 at which the date line crossing occurs.
*
* \returns the latitude at which the geodesic crosses the date line
* \since QGIS 3.6
*/
double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
double latitudeGeodesicCrossesDateLine( const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine SIP_OUT ) const;

private:

Expand Down
112 changes: 75 additions & 37 deletions tests/src/python/test_qgsdistancearea.py
Expand Up @@ -633,46 +633,84 @@ def testGeodesicIntersectionAtDateLine(self):
da.setSourceCrs(crs, QgsProject.instance().transformContext())
da.setEllipsoid("WGS84")

self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(0, 0), QgsPointXY(-170, 0)), 0, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-170, 0), QgsPointXY(170, 0)), 0, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(179, 0), QgsPointXY(181, 0)), 0, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-170, 0), QgsPointXY(170, 0)), 0, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(180, 0), QgsPointXY(180, 0)), 0, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(180, -10), QgsPointXY(180, -10)), -10, 5)

self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(138.26237, -20.314687), QgsPointXY(-151.6, -77.8)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(138.26237, -20.314687), QgsPointXY(-151.6 + 360, -77.8)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-151.6, -77.8), QgsPointXY(138.26237, -20.314687)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(170.60188754234980024, -70.81368329001529105),
QgsPointXY(-164.61259948055175073, -76.66761193248410677)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-164.61259948055175073, -76.66761193248410677),
QgsPointXY(170.60188754234980024,
-70.81368329001529105)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(178.44469761238570982, -73.47820480021761114),
QgsPointXY(-179.21026002627399976, -74.08952948682963324)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-179.21026002627399976, -74.08952948682963324),
QgsPointXY(178.44469761238570982,
-73.47820480021761114)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(179.83103440731269984, -73.8481044794813215),
QgsPointXY(-179.93191793815378787, -73.90885909527753483)), -73.89148222666744914, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-179.93191793815378787, -73.90885909527753483),
QgsPointXY(179.83103440731269984, -73.8481044794813215)), -73.89148222666744914, 5)

self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(179.92498611649580198, 7.24703528617311754),
QgsPointXY(-178.20070563806575592, 16.09649962419504732)), 7.6112109902580265, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-178.20070563806575592, 16.09649962419504732),
QgsPointXY(179.92498611649580198, 7.24703528617311754)), 7.6112109902580265, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(175.76717768974583578, 8.93749416467257873),
QgsPointXY(-175.15030911497356669, 8.59851183021221033)), 8.80683758146703966, 5)
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-175.15030911497356669, 8.59851183021221033),
QgsPointXY(175.76717768974583578,
8.93749416467257873)), 8.80683758146703966, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(0, 0), QgsPointXY(-170, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-170, 0), QgsPointXY(170, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0.5, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(179, 0), QgsPointXY(181, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0.5, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-170, 0), QgsPointXY(170, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0.5, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(180, 0), QgsPointXY(180, 0))
self.assertAlmostEqual(lat, 0, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(180, -10), QgsPointXY(180, -10))
self.assertAlmostEqual(lat, -10, 5)

lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(138.26237, -20.314687), QgsPointXY(-151.6, -77.8))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.007113545719515548, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(138.26237, -20.314687), QgsPointXY(-151.6 + 360, -77.8))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.007113545719515548, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-151.6, -77.8), QgsPointXY(138.26237, -20.314687))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.007113545719515548, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(170.60188754234980024, -70.81368329001529105),
QgsPointXY(-164.61259948055175073, -76.66761193248410677))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.0879577697523441, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-164.61259948055175073, -76.66761193248410677),
QgsPointXY(170.60188754234980024,
-70.81368329001529105))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.0879577697523441, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(178.44469761238570982, -73.47820480021761114),
QgsPointXY(-179.21026002627399976, -74.08952948682963324))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.6713541474159178, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-179.21026002627399976, -74.08952948682963324),
QgsPointXY(178.44469761238570982,
-73.47820480021761114))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.6713541474159178, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(179.83103440731269984, -73.8481044794813215),
QgsPointXY(-179.93191793815378787, -73.90885909527753483))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.7135414998986486, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-179.93191793815378787, -73.90885909527753483),
QgsPointXY(179.83103440731269984, -73.8481044794813215))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.7135414998986486, 5)

lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(179.92498611649580198, 7.24703528617311754),
QgsPointXY(-178.20070563806575592, 16.09649962419504732))
self.assertAlmostEqual(lat, 7.6112109902580265, 5)
self.assertAlmostEqual(fract, 0.04111771567489498, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-178.20070563806575592, 16.09649962419504732),
QgsPointXY(179.92498611649580198, 7.24703528617311754))
self.assertAlmostEqual(lat, 7.6112109902580265, 5)
self.assertAlmostEqual(fract, 0.04111771567489498, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(175.76717768974583578, 8.93749416467257873),
QgsPointXY(-175.15030911497356669, 8.59851183021221033))
self.assertAlmostEqual(lat, 8.80683758146703966, 5)
self.assertAlmostEqual(fract, 0.46581637044475815, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-175.15030911497356669, 8.59851183021221033),
QgsPointXY(175.76717768974583578,
8.93749416467257873))
self.assertAlmostEqual(lat, 8.80683758146703966, 5)
self.assertAlmostEqual(fract, 0.46581637044475815, 5)

# calculation should be ellipsoid dependent!
da.setEllipsoid("Phobos2000")
self.assertAlmostEqual(da.latitudeGeodesicCrossesDateLine(QgsPointXY(-175.15030911497356669, 8.59851183021221033),
QgsPointXY(175.76717768974583578,
8.93749416467257873)), 8.836479503936307, 5)
lat, fract = da.latitudeGeodesicCrossesDateLine(QgsPointXY(-175.15030911497356669, 8.59851183021221033),
QgsPointXY(175.76717768974583578,
8.93749416467257873))
self.assertAlmostEqual(lat, 8.836479503936307, 5)
self.assertAlmostEqual(fract, 0.46593364650414865, 5)

def testGeodesicLine(self):
da = QgsDistanceArea()
Expand Down

0 comments on commit 13ac907

Please sign in to comment.