Skip to content

Commit

Permalink
Fix some corner cases
Browse files Browse the repository at this point in the history
  • Loading branch information
nyalldawson committed Jan 11, 2019
1 parent 2a774d6 commit 629de0c
Show file tree
Hide file tree
Showing 5 changed files with 42 additions and 8 deletions.
Binary file not shown.
Binary file modified python/plugins/processing/tests/testdata/custom/pol.gpkg
Binary file not shown.
Binary file not shown.
10 changes: 10 additions & 0 deletions src/core/qgsdistancearea.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -476,6 +476,14 @@ double QgsDistanceArea::latitudeGeodesicCrossesAntimeridian( const QgsPointXY &p
double lat = p2y;
double lon = p2x;

if ( mEllipsoid == GEO_NONE )
{
fractionAlongLine = ( 180 - p1x ) / ( p2x - p1x );
if ( p1.x() >= 180 )
fractionAlongLine = 1 - fractionAlongLine;
return p1y + ( 180 - p1x ) / ( p2x - p1x ) * ( p2y - p1y );
}

geod_geodesic geod;
geod_init( &geod, mSemiMajor, 1 / mInvFlattening );

Expand Down Expand Up @@ -528,6 +536,8 @@ double QgsDistanceArea::latitudeGeodesicCrossesAntimeridian( const QgsPointXY &p
}

fractionAlongLine = intersectionDist / totalDist;
if ( p1.x() >= 180 )
fractionAlongLine = 1 - fractionAlongLine;

// either converged on 180 longitude or hit too many iterations
return lat;
Expand Down
40 changes: 32 additions & 8 deletions tests/src/python/test_qgsdistancearea.py
Original file line number Diff line number Diff line change
Expand Up @@ -649,6 +649,12 @@ def testGeodesicIntersectionAtAntimeridian(self):
self.assertAlmostEqual(lat, 0, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(180, -10), QgsPointXY(180, -10))
self.assertAlmostEqual(lat, -10, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(171, 0), QgsPointXY(181, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0.9, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(181, 0), QgsPointXY(171, 0))
self.assertAlmostEqual(lat, 0, 5)
self.assertAlmostEqual(fract, 0.1, 5)

lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(138.26237, -20.314687), QgsPointXY(-151.6, -77.8))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
Expand All @@ -658,7 +664,7 @@ def testGeodesicIntersectionAtAntimeridian(self):
self.assertAlmostEqual(fract, 0.007113545719515548, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(-151.6, -77.8), QgsPointXY(138.26237, -20.314687))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.007113545719515548, 5)
self.assertAlmostEqual(fract, 0.9928864542804845, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(170.60188754234980024, -70.81368329001529105),
QgsPointXY(-164.61259948055175073, -76.66761193248410677))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
Expand All @@ -667,7 +673,7 @@ def testGeodesicIntersectionAtAntimeridian(self):
QgsPointXY(170.60188754234980024,
-70.81368329001529105))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.0879577697523441, 5)
self.assertAlmostEqual(fract, 0.9120422302476558, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(178.44469761238570982, -73.47820480021761114),
QgsPointXY(-179.21026002627399976, -74.08952948682963324))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
Expand All @@ -676,15 +682,15 @@ def testGeodesicIntersectionAtAntimeridian(self):
QgsPointXY(178.44469761238570982,
-73.47820480021761114))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.6713541474159178, 5)
self.assertAlmostEqual(fract, 0.3286458525840822, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(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.latitudeGeodesicCrossesAntimeridian(QgsPointXY(-179.93191793815378787, -73.90885909527753483),
QgsPointXY(179.83103440731269984, -73.8481044794813215))
self.assertAlmostEqual(lat, -73.89148222666744914, 5)
self.assertAlmostEqual(fract, 0.7135414998986486, 5)
self.assertAlmostEqual(fract, 0.28645850010135143, 5)

lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(179.92498611649580198, 7.24703528617311754),
QgsPointXY(-178.20070563806575592, 16.09649962419504732))
Expand All @@ -693,11 +699,11 @@ def testGeodesicIntersectionAtAntimeridian(self):
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(-178.20070563806575592, 16.09649962419504732),
QgsPointXY(179.92498611649580198, 7.24703528617311754))
self.assertAlmostEqual(lat, 7.6112109902580265, 5)
self.assertAlmostEqual(fract, 0.04111771567489498, 5)
self.assertAlmostEqual(fract, 0.958882284325105, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(360 - 178.20070563806575592, 16.09649962419504732),
QgsPointXY(179.92498611649580198, 7.24703528617311754))
self.assertAlmostEqual(lat, 7.6112109902580265, 5)
self.assertAlmostEqual(fract, 0.04111771567489498, 5)
self.assertAlmostEqual(fract, 0.95888228432510, 5)

lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(175.76717768974583578, 8.93749416467257873),
QgsPointXY(-175.15030911497356669, 8.59851183021221033))
Expand All @@ -707,15 +713,33 @@ def testGeodesicIntersectionAtAntimeridian(self):
QgsPointXY(175.76717768974583578,
8.93749416467257873))
self.assertAlmostEqual(lat, 8.80683758146703966, 5)
self.assertAlmostEqual(fract, 0.46581637044475815, 5)
self.assertAlmostEqual(fract, 0.5341836295552418, 5)

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

# no ellipsoid
da.setEllipsoid("NONE")
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(-175, 8),
QgsPointXY(175,
9))
self.assertAlmostEqual(lat, 8.5, 5)
self.assertAlmostEqual(fract, 0.5, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(165, 8),
QgsPointXY(-175,
9))
self.assertAlmostEqual(lat, 8.75, 5)
self.assertAlmostEqual(fract, 0.75, 5)
lat, fract = da.latitudeGeodesicCrossesAntimeridian(QgsPointXY(-175, 8),
QgsPointXY(165,
9))
self.assertAlmostEqual(lat, 8.25, 5)
self.assertAlmostEqual(fract, 0.25, 5)

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

0 comments on commit 629de0c

Please sign in to comment.