Fix some corner cases

qgis · Jan 11, 2019 · 629de0c · 629de0c
1 parent 2a774d6
commit 629de0c
Show file tree

Hide file tree

Showing 5 changed files with 42 additions and 8 deletions.
diff --git a/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg b/python/plugins/processing/tests/testdata/custom/circular_strings.gpkg
diff --git a/python/plugins/processing/tests/testdata/custom/pol.gpkg b/python/plugins/processing/tests/testdata/custom/pol.gpkg
diff --git a/python/plugins/processing/tests/testdata/expected/merged_pol.gpkg b/python/plugins/processing/tests/testdata/expected/merged_pol.gpkg
diff --git a/src/core/qgsdistancearea.cpp b/src/core/qgsdistancearea.cpp
@@ -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 );
 
@@ -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;

diff --git a/tests/src/python/test_qgsdistancearea.py b/tests/src/python/test_qgsdistancearea.py
@@ -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)
@@ -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)
@@ -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)
@@ -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))
@@ -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))
@@ -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()