libgis: improve and explain fix for G_ellipsoid_polygon_area()

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@71169 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/gis/area_poly1.c

@@ -131,9 +131,7 @@ double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
     double x1, y1, x2, y2, dx, dy;
     double Qbar1, Qbar2;
     double area;
-    double thresh = 0.1;	/* threshold for latitude differences in arc seconds */
-    
-    thresh = Radians(thresh / 3600.0);
+    double thresh = 1e-6;	/* threshold for dy, should be between 1e-4 and 1e-7 */
 
     x2 = Radians(lon[n - 1]);
     y2 = Radians(lat[n - 1]);
@@ -163,9 +161,18 @@ double G_ellipsoid_polygon_area(const double *lon, const double *lat, int n)
 	if (fabs(dy) > thresh) {
 	    /* account for different latitudes y1, y2 */
 	    area += dx * (st->Qp - (Qbar2 - Qbar1) / dy);
+	    /* original: 
+	     * area += dx * st->Qp - (dx / dy) * (Qbar2 - Qbar1);
+	     */
 	}
 	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 * (st->Qp - Q(y2));
 	}
     }