diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index 3fc59eefe8f..cb9a24c341e 100644 --- a/lib/gis/area_poly2.c +++ b/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); } diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index 74e1622cc70..ddd47a4e7e4 100644 --- a/lib/vector/diglib/poly.c +++ b/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; diff --git a/vector/v.build/testsuite/data/tiny_polygons.pack b/vector/v.build/testsuite/data/tiny_polygons.pack new file mode 100644 index 00000000000..167516a8c7b Binary files /dev/null and b/vector/v.build/testsuite/data/tiny_polygons.pack differ diff --git a/vector/v.build/testsuite/test_v_build.py b/vector/v.build/testsuite/test_v_build.py index 03915edba86..3f322ef7cea 100644 --- a/vector/v.build/testsuite/test_v_build.py +++ b/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() diff --git a/vector/v.to.db/testsuite/test_v_to_db.py b/vector/v.to.db/testsuite/test_v_to_db.py index 723af416de7..8516eaaa350 100644 --- a/vector/v.to.db/testsuite/test_v_to_db.py +++ b/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])