v.overlay: use Vect_snap_line()

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@54132 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/v.overlay/area_area.c

@@ -15,6 +15,15 @@
 #include <grass/glocale.h>
 #include "local.h"
 
+/* for ilist qsort'ing and bsearch'ing */
+static int cmp_int(const void *a, const void *b)
+{
+    int ai = *(int *)a;
+    int bi = *(int *)b;
+    
+    return (ai < bi ? -1 : (ai > bi));
+}
+
 int area_area(struct Map_info *In, int *field, struct Map_info *Out,
 	      struct field_info *Fi, dbDriver * driver, int operator,
 	      int *ofield, ATTRIBUTES * attr, struct ilist *BList, double snap)
@@ -34,13 +43,63 @@ int area_area(struct Map_info *In, int *field, struct Map_info *Out,
 
     Points = Vect_new_line_struct();
     Cats = Vect_new_cats_struct();
-    
+
     /* optional snap */
     if (snap > 0) {
-	G_message(_("Snapping lines..."));
+	int i, j, snapped_lines = 0;
+	struct bound_box box;
+	struct boxlist *boxlist = Vect_new_boxlist(0);
+	struct ilist *reflist = Vect_new_list();
+	
+	G_message(_("Snapping boundaries with %g ..."), snap);
+
 	/* snap boundaries in B to boundaries in A
 	 * not modifying boundaries in A */
-	Vect_snap_lines_list2(Out, BList, GV_BOUNDARY, snap, NULL);
+
+	if (BList->n_values > 1)
+	    qsort(BList->value, BList->n_values, sizeof(int), cmp_int);
+
+	snapped_lines = 0;
+	nlines = BList->n_values;
+	for (i = 0; i < nlines; i++) {
+	    line = BList->value[i];
+	    Vect_read_line(Out, Points, Cats, line);
+	    /* select lines by box */
+	    Vect_get_line_box(Out, line, &box);
+	    box.E += snap;
+	    box.W -= snap;
+	    box.N += snap;
+	    box.S -= snap;
+	    box.T = 0.0;
+	    box.B = 0.0;
+	    Vect_select_lines_by_box(Out, &box, GV_BOUNDARY, boxlist);
+	    
+	    if (boxlist->n_values > 0) {
+		Vect_reset_list(reflist);
+		for (j = 0; j < boxlist->n_values; j++) {
+		    int aline = boxlist->id[j];
+
+		    if (!bsearch(&aline, BList->value, BList->n_values,
+			sizeof(int), cmp_int)) {
+			G_ilist_add(reflist, aline);
+		    }
+		}
+		
+		/* snap bline to alines */
+		if (Vect_snap_line(Out, reflist, Points, snap, 0, NULL, NULL)) {
+		    /* rewrite bline*/
+		    Vect_delete_line(Out, line);
+		    ret = Vect_write_line(Out, GV_BOUNDARY, Points, Cats);
+		    G_ilist_add(BList, ret);
+		    snapped_lines++;
+		    G_debug(3, "line %d snapped", line);
+		}
+	    }
+	}
+	Vect_destroy_boxlist(boxlist);
+	Vect_destroy_list(reflist);
+
+	G_verbose_message(_("%d boundaries snapped"), snapped_lines);
     }
 
     /* same procedure like for v.in.ogr
@@ -49,7 +108,7 @@ int area_area(struct Map_info *In, int *field, struct Map_info *Out,
      * and Vect_clean_small_angles_at_nodes() until no more small dangles are found */
     do {
 	G_message(_("Breaking lines..."));
-	Vect_break_lines_list(Out, NULL, BList, GV_LINE | GV_BOUNDARY, NULL);
+	Vect_break_lines_list(Out, NULL, BList, GV_BOUNDARY, NULL);
 
 	/* Probably not necessary for LINE x AREA */
 	G_message(_("Removing duplicates..."));

vector/v.overlay/main.c

@@ -9,6 +9,7 @@
  *               Markus Neteler <neteler itc.it>,
  *               Paul Kelly <paul-grass stjohnspoint.co.uk>
  *               OGR support by Martin Landa <landa.martin gmail.com>
+ *               Markus Metz
  * PURPOSE:      
  * COPYRIGHT:    (C) 2003-2009 by the GRASS Development Team
  *
@@ -228,6 +229,7 @@ int main(int argc, char *argv[])
 	nlines_out = 0;
 	for (line = 1; line <= nlines; line++) {
 	    int ltype;
+	    int vertices = 100; /* max number of vertices per line */
 
 	    G_percent(line, nlines, 1);	/* must be before any continue */
 
@@ -248,10 +250,11 @@ int main(int argc, char *argv[])
 		continue;
 
 	    /* TODO: figure out a reasonable threshold */
-	    if (Points->n_points > 100) {
-		int vertices = 100;
+	    if (Points->n_points > vertices) {
 		int start = 0;	/* number of coordinates written */
 		
+		vertices = Points->n_points / (Points->n_points / vertices + 1);
+		
 		/* split */
 		while (start < Points->n_points - 1) {
 		    int v = 0;

vector/v.overlay/v.overlay.html

@@ -11,9 +11,12 @@ it is sometimes necessary to snap areas of <em>binput</em> to those of
 can be created. Snapping is enabled by default and can be disabled by 
 setting the <em>snap</em> option to a negative value. Recommended values 
 are between 0.00000001 and 0.0001. Using larger values for snapping can 
-have undesired side-effects. Snapping modifies only boundaries in binput, 
-which are snapped to boundaries in ainput. Boundaries in ainput are not 
-modified.
+have undesired side-effects, but may sometimes be necessary to get a 
+clean output (see example below). In general, it is recommended to start 
+with a small snapping threshold, gradually increasing the threshold until 
+the result is reasonably clean. Snapping modifies only boundaries in 
+binput, which are snapped to boundaries in ainput. Boundaries in ainput 
+are not modified.
 <!-- This is outdated 
 There are 3 links attached to features in output map, 
 <ul>
@@ -70,13 +73,13 @@ v.db.connect map=outputmap table=binput.dbf field=3
 it! Therefore it is advisable to copy tables from ainput and binput first and
 connect the copied tables to the output map.-->
 
-<h2>EXAMPLE POLYGON TO POLYGON UNION</h2>
+<h2>EXAMPLES</h2>
+
+<h4>Polygons overlaid with ploygons</h4>
 <div class="code"><pre>
 v.overlay ainput=lake binput=province output=lakeXprovince operator=or
 </pre></div>
 
-<h2>EXAMPLE POLYGON TO POLYGON UNION</h2>
-
 Polygon union of urban area and Census 2000 areas (North Carolina dataset):
 
 <div class="code"><pre>
@@ -120,11 +123,14 @@ b_UA_TYPE|UA
 <img src="v_overlay_census_wake2000.png" alt="GRASS v.overlay: polygon to polygon union (input 2)" border=1>
 <img src="v_overlay_urban_census2000.png" alt="GRASS v.overlay: polygon to polygon union (result)" border=1>
 <br>
-<i>v.overlay: Polygon union of urban area and Census 2000 areas (North Carolina dataset)</i>
+<i>v.overlay: Polygon union (right) of urban area (left) and Census 2000 (middle) areas (North Carolina dataset)</i>
 </center>
 
+As can be seen by the resulting large number of centroids on boundaries, 
+the urban areas do not match exactly the Census 2000 areas. In this case 
+a clean result can be obtained by snapping with a threshold of 0.1 (1 cm).
 
-<h2>EXAMPLE LINE TO POLYGON CLIPPING</h2>
+<h4>Lines overlaid with polygons</h4>
 
 Using the North Carolina sample dataset, we clip the roads map to the area
 of city of Raleigh, preserving road attributes in layer 1: