From d85039363a552b901da81467f45925f183bd50a6 Mon Sep 17 00:00:00 2001
From: Nyall Dawson <nyall.dawson@gmail.com>
Date: Mon, 10 Jul 2017 09:31:16 +1000
Subject: [PATCH] Fix incorrect area calculation in corner cases (fix #16820)

In certain circumstances very proximal nodes could cause instability
in the ellipsoidal area calculation.

Port (slightly tweaked) fix from grass changeset 71167 for same issue,
and add a unit test
---
 src/core/qgsdistancearea.cpp           | 33 ++++++++++++++++++++++----
 tests/src/core/testqgsdistancearea.cpp | 11 +++++++++
 2 files changed, 40 insertions(+), 4 deletions(-)

diff --git a/src/core/qgsdistancearea.cpp b/src/core/qgsdistancearea.cpp
index a9e4255af19..8450eba39db 100644
--- a/src/core/qgsdistancearea.cpp
+++ b/src/core/qgsdistancearea.cpp
@@ -678,6 +678,14 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
   double Qbar1, Qbar2;
   double area;
 
+  /* GRASS comment: threshold for dy, should be between 1e-4 and 1e-7
+   * QGIS note: while the grass comment states that thres should be between 1e-4->1e-7,
+   * a value of 1e-7 caused TestQgsDistanceArea::regression14675() to regress
+   * The maximum threshold possible which permits regression14675() to pass
+   * was found to be ~0.7e-7.
+  */
+  const double thresh = 0.7e-7;
+
   QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
   if ( !willUseEllipsoid() )
   {
@@ -708,11 +716,28 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
         x1 += m_TwoPI;
 
     dx = x2 - x1;
-    area += dx * ( m_Qp - getQ( y2 ) );
-
     dy = y2 - y1;
-    if ( !qgsDoubleNear( dy, 0.0 ) )
-      area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
+    if ( qAbs( dy ) > thresh )
+    {
+      /* account for different latitudes y1, y2 */
+      area += dx * ( m_Qp - ( Qbar2 - Qbar1 ) / dy );
+    }
+    else
+    {
+      /* latitudes y1, y2 are (nearly) identical */
+
+      /* if y2 becomes similar to y1, i.e. y2 -> y1
+       * Qbar2 - Qbar1 -> 0 and dy -> 0
+       * (Qbar2 - Qbar1) / dy -> ?
+       * (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
+       * Metz 2017
+       */
+      area += dx * ( m_Qp - getQ( y2 ) );
+
+      /* original:
+       * area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
+       */
+    }
   }
   if ( ( area *= m_AE ) < 0.0 )
     area = -area;
diff --git a/tests/src/core/testqgsdistancearea.cpp b/tests/src/core/testqgsdistancearea.cpp
index 463752d9936..4286c64e3cd 100644
--- a/tests/src/core/testqgsdistancearea.cpp
+++ b/tests/src/core/testqgsdistancearea.cpp
@@ -46,6 +46,7 @@ class TestQgsDistanceArea: public QObject
     void measureAreaAndUnits();
     void emptyPolygon();
     void regression14675();
+    void regression16820();
 
 };
 
@@ -379,6 +380,16 @@ void TestQgsDistanceArea::regression14675()
   QGSCOMPARENEAR( calc.measureArea( geom ), 0.83301, 0.02 );
 }
 
+void TestQgsDistanceArea::regression16820()
+{
+  QgsDistanceArea calc;
+  calc.setEllipsoid( QStringLiteral( "WGS84" ) );
+  calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:32634" ) ) );
+  QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((110250.54038314701756462 5084495.57398066483438015, 110243.46975068224128336 5084507.17200060561299324, 110251.23908144699817058 5084506.68309532757848501, 110251.2394439501222223 5084506.68307251576334238, 110250.54048078990308568 5084495.57553235255181789, 110250.54038314701756462 5084495.57398066483438015))" ) ) );
+  //lots of tolerance here - the formulas get quite unstable with small areas due to division by very small floats
+  QGSCOMPARENEAR( calc.measureArea( geom ), 43.183369, 0.02 );
+}
+
 QGSTEST_MAIN( TestQgsDistanceArea )
 #include "testqgsdistancearea.moc"