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

[svn_grass at osgeo.org]

Author: mmetz
Date: 2017-06-07 05:59:08 -0700 (Wed, 07 Jun 2017)
New Revision: 71167

Modified:
   grass/trunk/lib/gis/area_poly1.c
Log:
libgis: handle subtle differences in G_ellipsoid_polygon_area() (#3356)

Modified: grass/trunk/lib/gis/area_poly1.c
===================================================================
--- grass/trunk/lib/gis/area_poly1.c	2017-06-05 18:02:55 UTC (rev 71166)
+++ grass/trunk/lib/gis/area_poly1.c	2017-06-07 12:59:08 UTC (rev 71167)
@@ -131,6 +131,9 @@
     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);
 
     x2 = Radians(lon[n - 1]);
     y2 = Radians(lat[n - 1]);
@@ -155,10 +158,16 @@
 		x1 += TWOPI;
 
 	dx = x2 - x1;
-	area += dx * (st->Qp - Q(y2));
 
-	if ((dy = y2 - y1) != 0.0)
-	    area += dx * Q(y2) - (dx / dy) * (Qbar2 - Qbar1);
+	dy = y2 - y1;
+	if (fabs(dy) > thresh) {
+	    /* account for different latitudes y1, y2 */
+	    area += dx * (st->Qp - (Qbar2 - Qbar1) / dy);
+	}
+	else {
+	    /* latitudes y1, y2 are (nearly) identical */
+	    area += dx * (st->Qp - Q(y2));
+	}
     }
     if ((area *= st->AE) < 0.0)
 	area = -area;