[GRASS-SVN] r71169 - grass/trunk/lib/gis

[svn_grass at osgeo.org]

Author: mmetz
Date: 2017-06-07 13:14:16 -0700 (Wed, 07 Jun 2017)
New Revision: 71169

Modified:
   grass/trunk/lib/gis/area_poly1.c
Log:
libgis: improve and explain fix for G_ellipsoid_polygon_area()

Modified: grass/trunk/lib/gis/area_poly1.c
===================================================================
--- grass/trunk/lib/gis/area_poly1.c	2017-06-07 13:38:25 UTC (rev 71168)
+++ grass/trunk/lib/gis/area_poly1.c	2017-06-07 20:14:16 UTC (rev 71169)
@@ -131,9 +131,7 @@
     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 @@
 	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));
 	}
     }