v.overlay: add option to remove small areas

vector/v.overlay/area_area.c

@@ -26,7 +26,7 @@ static int cmp_int(const void *a, const void *b)
 int area_area(struct Map_info *In, int *field, struct Map_info *Tmp,
               struct Map_info *Out, struct field_info *Fi, dbDriver *driver,
               int operator, int *ofield, ATTRIBUTES *attr, struct ilist *BList,
-              double snap)
+              double snap, double area_minsize)
 {
     int ret, input, line, nlines, area, nareas;
     int in_centr, out_cat;
@@ -172,10 +172,19 @@ int area_area(struct Map_info *In, int *field, struct Map_info *Tmp,
         Vect_remove_bridges(Tmp, NULL, NULL, NULL);
     }
 
-    G_set_verbose(0);
-    Vect_build_partial(Tmp, GV_BUILD_NONE);
-    Vect_build_partial(Tmp, GV_BUILD_BASE);
-    G_set_verbose(verbose);
+    if (area_minsize > 0) {
+        Vect_build_partial(Tmp, GV_BUILD_CENTROIDS);
+        G_message(_("Removing areas smaller than %g sqm..."), area_minsize);
+        Vect_remove_small_areas(Tmp, area_minsize, NULL, NULL);
+        Vect_build_partial(Tmp, GV_BUILD_BASE);
+    }
+    else {
+        G_set_verbose(0);
+        Vect_build_partial(Tmp, GV_BUILD_NONE);
+        Vect_build_partial(Tmp, GV_BUILD_BASE);
+        G_set_verbose(verbose);
+    }
+
     G_message(_("Merging lines..."));
     Vect_merge_lines(Tmp, GV_BOUNDARY, NULL, NULL);
 

vector/v.overlay/local.h

@@ -29,7 +29,7 @@ ATTR *find_attr(ATTRIBUTES *attributes, int cat);
 int area_area(struct Map_info *In, int *field, struct Map_info *Tmp,
               struct Map_info *Out, struct field_info *Fi, dbDriver *driver,
               int operator, int * ofield, ATTRIBUTES *attr, struct ilist *BList,
-              double snap_thresh);
+              double snap_thresh, double area_minsize);
 int line_area(struct Map_info *In, int *field, struct Map_info *Tmp,
               struct Map_info *Out, struct field_info *Fi, dbDriver *driver,
               int operator, int * ofield, ATTRIBUTES *attr,

vector/v.overlay/main.c

@@ -10,7 +10,7 @@
  *               OGR support by Martin Landa <landa.martin gmail.com>
  *               Markus Metz
  * PURPOSE:
- * COPYRIGHT:    (C) 2003-2016 by the GRASS Development Team
+ * COPYRIGHT:    (C) 2003-2026 by the GRASS Development Team
  *
  *               This program is free software under the GNU General
  *               Public License (>=v2). Read the file COPYING that
@@ -36,10 +36,10 @@ int main(int argc, char *argv[])
 {
     int i, j, input, line, nlines, operator;
     int type[2], field[2], ofield[3];
-    double snap_thresh;
+    double snap_thresh, area_minsize;
     struct GModule *module;
     struct Option *in_opt[2], *out_opt, *type_opt[2], *field_opt[2],
-        *ofield_opt, *operator_opt, *snap_opt;
+        *ofield_opt, *operator_opt, *snap_opt, *minsize_opt;
     struct Flag *table_flag;
     struct Map_info In[2], Out, Tmp;
     struct line_pnts *Points, *Points2;
@@ -142,6 +142,13 @@ int main(int argc, char *argv[])
     snap_opt->type = TYPE_DOUBLE;
     snap_opt->answer = "1e-8";
 
+    minsize_opt = G_define_option();
+    minsize_opt->key = "minsize";
+    minsize_opt->label = _("Minimum size of areas to keep in square meters");
+    minsize_opt->description = _("Enable with minsize > 0");
+    minsize_opt->type = TYPE_DOUBLE;
+    minsize_opt->answer = "0";
+
     table_flag = G_define_standard_flag(G_FLG_V_TABLE);
     table_flag->guisection = _("Attributes");
 
@@ -220,6 +227,7 @@ int main(int argc, char *argv[])
                       operator_opt->answer);
 
     snap_thresh = atof(snap_opt->answer);
+    area_minsize = atof(minsize_opt->answer);
 
     Points = Vect_new_line_struct();
     Points2 = Vect_new_line_struct();
@@ -635,7 +643,7 @@ int main(int argc, char *argv[])
     /* AREA x AREA */
     if (type[0] == GV_AREA) {
         area_area(In, field, &Tmp, &Out, Fi, driver, operator, ofield, attr,
-                  BList, snap_thresh);
+                  BList, snap_thresh, area_minsize);
     }
     else { /* LINE x AREA */
         line_area(In, field, &Tmp, &Out, Fi, driver, operator, ofield, attr,

vector/v.overlay/tests/test_v_overlay.py

@@ -0,0 +1,128 @@
+import pytest
+import grass.script as gs
+
+# run in a GRASS session with nc_spm_full_v2beta1
+
+
+class TestVOverlay:
+    """Test v.overlay output against expected output"""
+
+    # create test data
+    @pytest.fixture(scope="class", autouse=True)
+    def create_testdata(self):
+        # set up
+        gs.run_command(
+            "v.extract",
+            input="boundary_county",
+            output="boundary_county_extract1",
+            where="NAME in ('CURRITUCK')",
+        )
+        gs.run_command(
+            "v.extract",
+            input="boundary_county",
+            output="boundary_county_extract2",
+            where="NAME in ('CAMDEN')",
+        )
+        # modify extract 1
+        gs.run_command(
+            "v.buffer",
+            input="boundary_county_extract1",
+            output="boundary_county_extract1_buffer_out",
+            type="area",
+            distance=2,
+        )
+        gs.run_command(
+            "v.buffer",
+            input="boundary_county_extract1_buffer_out",
+            output="boundary_county_extract1_buffer_in",
+            type="area",
+            distance=-2,
+        )
+
+        # run the tests
+        yield
+
+        # clean up test data regardless of test success/failure
+        gs.run_command(
+            "g.remove", type="vector", flags="f", pattern="boundary_county_extract*"
+        )
+
+    def test_voverlay_nocleaning(self):
+        """
+        Overlay two input vectors and compare the output to the expected output.
+        The output will have many very small areas.
+        This test would fail in GRASS84 with a too small number of centroids.
+        Since GRASS84, calculation for centroids has been improved for
+        edge cases like very small areas."""
+        gs.run_command(
+            "v.overlay",
+            ainput="boundary_county_extract1_buffer_in",
+            atype="area",
+            binput="boundary_county_extract2",
+            btype="area",
+            output="boundary_county_extract_overlay_nocleaning",
+            operator="or",
+            snap=-1,
+            quiet=True,
+        )
+
+        reference = {
+            "nodes": 1050,
+            "primitives": 1802,
+            "points": 0,
+            "lines": 0,
+            "boundaries": 1572,
+            "centroids": 230,
+            "areas": 529,
+            "islands": 7,
+        }
+
+        observed = gs.vector_info_topo(
+            map="boundary_county_extract_overlay_nocleaning",
+        )
+
+        # compare results
+        for key in list(reference):
+            # all integers, simple comparison is ok
+            assert reference[key] == observed[key], (
+                f"Difference in {key}, expected {reference[key]}, got {observed[key]}"
+            )
+
+    def test_voverlay_withcleaning(self):
+        """Overlay two input vectors and compare the output to the expected output.
+        The output will have many very small areas.
+        This is a test for the cleaning options snap and minsize"""
+        gs.run_command(
+            "v.overlay",
+            ainput="boundary_county_extract1_buffer_in",
+            atype="area",
+            binput="boundary_county_extract2",
+            btype="area",
+            output="boundary_county_extract_overlay_withcleaning",
+            operator="or",
+            snap=1.0e-6,
+            minsize=0.00001,
+            quiet=True,
+        )
+
+        reference = {
+            "nodes": 1232,
+            "primitives": 2164,
+            "points": 0,
+            "lines": 0,
+            "boundaries": 1847,
+            "centroids": 317,
+            "areas": 622,
+            "islands": 7,
+        }
+
+        observed = gs.vector_info_topo(
+            map="boundary_county_extract_overlay_withcleaning",
+        )
+
+        # compare results
+        for key in list(reference):
+            # all integers, simple comparison is ok
+            assert reference[key] == observed[key], (
+                f"Difference in {key}, expected {reference[key]}, got {observed[key]}"
+            )

vector/v.overlay/v.overlay.html

@@ -55,8 +55,16 @@ <h2>NOTES</h2>
 the output layer with the third number.
 
 <p>
-If <b>atype</b>=auto is given than <em>v.overlay</em> determines
+If <b>atype</b>=auto is given then <em>v.overlay</em> determines the
 feature type for <b>ainput</b> from the first found feature.
+<p>
+When overlaying two vectors with areas, very small areas can occur in the
+output. This can happen when e.g. one vector is a slightly modified
+version of the other vector (buffered or simplified). These very small
+areas can be removed by setting <b>minsize</b> to some value larger 0.
+The value is interpreted as square meters. In order to remove only noise
+from slightly mismatching boundaries, the value of <b>minsize</b> should be
+small, e.g. in the range 0.0001 to 1.
 
 <!-- This is outdated
 <p><div class="code"><pre>

vector/v.overlay/v.overlay.md

@@ -37,8 +37,16 @@ If the *third* number of the **olayer** option is greater than 0, then
 the categories of **binput** in layer *blayer* are transferred to the
 output layer with the third number.
 
-If **atype**=auto is given than *v.overlay* determines feature type for
-**ainput** from the first found feature.
+If **atype**=auto is given then *v.overlay* determines the feature type
+for **ainput** from the first found feature.
+
+When overlaying two vectors with areas, very small areas can occur in the
+output. This can happen when e.g. one vector is a slightly modified
+version of the other vector (buffered or simplified). These very small
+areas can be removed by setting **minsize** to some value larger 0.
+The value is interpreted as square meters. In order to remove only noise
+from slightly mismatching boundaries, the value of **minsize** should be
+small, e.g. in the range 0.0001 to 1.
 
 ## EXAMPLES