lib: more robust signed polygon size

lib/gis/area_poly2.c

@@ -24,11 +24,31 @@
  */
 double G_planimetric_polygon_area(const double *x, const double *y, int n)
 {
-    double x1, y1, x2, y2;
+    double x1, y1, x2, y2, xshift, yshift;
     double area;
 
-    x2 = x[n - 1];
-    y2 = y[n - 1];
+    /* Get a reference point close to the polygon as new origin.
+     * The first point would be good enough, particularly for small
+     * polygons. The average of the first and the mid-point is a fast
+     * approximation of a reference point with reduced distance to the
+     * ring's vertices, also for larger rings.
+     *
+     * Shift coordinates towards this reference point to make
+     * calculation of the signed area more robust by increasing the
+     * accuracy for given fp precision limits.
+     *
+     * Considering the basic formula
+     *    area += (x2 - x1) * (y2 + y1)
+     * the shift is in theory only needed for addition, not subtraction,
+     * but does no harm and is kept for symmetry treating x and y coords.
+     *
+     * Keep in sync with dig_find_area_poly() in lib/vector/diglib/poly.c */
+
+    xshift = (x[0] + x[n / 2]) / 2.;
+    yshift = (y[0] + y[n / 2]) / 2.;
+
+    x2 = x[n - 1] - xshift;
+    y2 = y[n - 1] - yshift;
 
     area = 0;
     while (--n >= 0) {
@@ -38,6 +58,9 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n)
         x2 = *x++;
         y2 = *y++;
 
+        x2 -= xshift;
+        y2 -= yshift;
+
         area += (y2 + y1) * (x2 - x1);
     }
 

lib/vector/diglib/poly.c

@@ -96,19 +96,41 @@ int dig_get_poly_points(int n_lines, struct line_pnts **LPoints,
  */
 int dig_find_area_poly(struct line_pnts *Points, double *totalarea)
 {
-    int i;
-    double *x, *y;
+    int i, mid;
+    double *x, *y, xshift, yshift;
     double tot_area;
 
     x = Points->x;
     y = Points->y;
 
+    /* Get a reference point close to the polygon as new origin.
+     * The first point would be good enough, particularly for small
+     * polygons. The average of the first and the mid-point is a fast
+     * approximation of a reference point with reduced distance to the
+     * ring's vertices, also for larger rings.
+     *
+     * Shift coordinates towards this reference point to make
+     * calculation of the signed area more robust by increasing the
+     * accuracy for given fp precision limits.
+     *
+     * Considering the basic formula
+     *    area += (x2 - x1) * (y2 + y1)
+     * the shift is in theory only needed for addition, not subtraction,
+     * but does no harm and is kept for symmetry treating x and y coords.
+     *
+     * Keep in sync with G_planimetric_polygon_area() in lib/gis/area_poly2.c */
+
+    mid = Points->n_points / 2;
+    xshift = (x[0] + x[mid]) / 2.;
+    yshift = (y[0] + y[mid]) / 2.;
+
     /* line integral: *Points do not need to be pruned */
     /* surveyor's formula is more common, but more prone to
      * fp precision limit errors, and *Points would need to be pruned */
     tot_area = 0.0;
     for (i = 1; i < Points->n_points; i++) {
-        tot_area += (x[i] - x[i - 1]) * (y[i] + y[i - 1]);
+        tot_area += ((x[i] - xshift) - (x[i - 1] - xshift)) *
+                    ((y[i] - yshift) + (y[i - 1] - yshift));
     }
     *totalarea = 0.5 * tot_area;
 

vector/v.build/testsuite/test_v_build.py

@@ -6,6 +6,8 @@
 class TestVBuild(TestCase):
     """Test v.build output against expected output stored in vbuild_output"""
 
+    tiny_polygons = "test_tiny_polygons"
+
     @classmethod
     def setUpClass(cls):
         """Set up a temporary region and create vector map with test features."""
@@ -64,10 +66,25 @@ def setUpClass(cls):
         0 |    2 |         2
 ------------------------------------------------------------------------------------------"""
 
+        cls.runModule(
+            "v.unpack",
+            input="data/tiny_polygons.pack",
+            output=cls.tiny_polygons,
+        )
+        # This test vector contains a few small areas. One of the areas
+        # is so tiny that the sign of the signed area size might be wrong
+        # due to fp precision limits. In this case, the area is
+        # erroneously identified as an island and not as an area
+        cls.runModule(
+            "v.build",
+            map=cls.tiny_polygons,
+        )
+
     @classmethod
     def tearDownClass(cls):
         """Remove created vector map and temporary files, then delete the temp region."""
         gs.run_command("g.remove", type="vector", flags="f", name="test_3x3_map")
+        gs.run_command("g.remove", type="vector", flags="f", name=cls.tiny_polygons)
         cls.del_temp_region()
 
     def test_vbuild_output(self):
@@ -87,6 +104,26 @@ def test_vbuild_output(self):
         )
         self.assertMultiLineEqual(filtered_output, self.vbuild_output)
 
+    def test_areas_islands(self):
+        """Compare the v.info output to the expected output."""
+        self.assertModuleKeyValue(
+            "v.info",
+            map=self.tiny_polygons,
+            flags="tg",
+            sep="=",
+            precision=0.1,
+            reference={
+                "nodes": 159,
+                "points": 0,
+                "lines": 0,
+                "boundaries": 168,
+                "centroids": 5,
+                "areas": 10,
+                "islands": 1,
+                "primitives": 173,
+            },
+        )
+
 
 if __name__ == "__main__":
     test()

vector/v.to.db/testsuite/test_v_to_db.py

@@ -1,4 +1,5 @@
 import json
+import math
 from itertools import zip_longest
 
 from grass.gunittest.case import TestCase
@@ -16,7 +17,17 @@ def _assert_dict_almost_equal(self, d1, d2):
         self.assertEqual(set(d1.keys()), set(d2.keys()))
         for k1 in d1:
             if isinstance(d1[k1], float):
-                self.assertAlmostEqual(d1[k1], d2[k1], places=6)
+                m_obs, e_obs = math.frexp(d1[k1])
+                m_ref, e_ref = math.frexp(d2[k1])
+                self.assertEqual(e_obs, e_ref)
+                # use 12 instead of default 7 decimal places
+                # for reliable double precision fp results
+                self.assertAlmostEqual(
+                    m_obs,
+                    m_ref,
+                    places=12,
+                    msg=f"Comparing mantissas of {d1[k1]} and {d2[k1]}",
+                )
             else:
                 self.assertEqual(d1[k1], d2[k1])