vectorintro corrected

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

vector/v.in.ogr/geom.c

@@ -25,6 +25,7 @@
 #include "ogr_api.h"
 #include "global.h"
 
+int split_line(struct Map_info *Map, int otype, struct line_pnts *Points, struct line_cats *Cats);
 
 /* Add categories to centroids inside polygon */
 int
@@ -206,7 +207,11 @@ geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 	    otype = GV_BOUNDARY;
 	else
 	    otype = GV_LINE;
-	Vect_write_line(Map, otype, Points, Cats);
+
+	if (split_distance > 0 && otype == GV_BOUNDARY)
+	    split_line(Map, otype, Points, Cats);
+	else
+	    Vect_write_line(Map, otype, Points, Cats);
     }
 
     else if (eType == wkbPolygon) {
@@ -243,7 +248,11 @@ geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 	    otype = GV_LINE;
 	else
 	    otype = GV_BOUNDARY;
-	Vect_write_line(Map, otype, Points, BCats);
+
+	if (split_distance > 0 && otype == GV_BOUNDARY)
+	    split_line(Map, otype, Points, Cats);
+	else
+	    Vect_write_line(Map, otype, Points, Cats);
 
 	/* Isles */
 	IPoints =
@@ -274,7 +283,11 @@ geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 		    otype = GV_LINE;
 		else
 		    otype = GV_BOUNDARY;
-		Vect_write_line(Map, otype, IPoints[i - 1], BCats);
+
+		if (split_distance > 0 && otype == GV_BOUNDARY)
+		    split_line(Map, otype, IPoints[i - 1], BCats);
+		else
+		    Vect_write_line(Map, otype, IPoints[i - 1], BCats);
 	    }
 	}
 
@@ -349,3 +362,34 @@ geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 
     return 0;
 }
+
+int split_line(struct Map_info *Map, int otype, struct line_pnts *Points, struct line_cats *Cats)
+{
+    int i;
+    double dist = 0., dx, dy;
+    struct line_pnts *OutPoints;
+
+    OutPoints = Vect_new_line_struct();
+    Vect_reset_line(OutPoints);
+    i = 0;
+    Vect_append_point(OutPoints, Points->x[i], Points->y[i], Points->z[i]);
+
+    for (i = 1; i < Points->n_points; i++) {
+	dx = Points->x[i] - Points->x[i - 1];
+	dy = Points->y[i] - Points->y[i - 1];
+	dist += sqrt(dx * dx + dy * dy);
+	if (OutPoints->n_points > 1 && dist > split_distance) {
+	    Vect_write_line(Map, otype, OutPoints, Cats);
+	    Vect_reset_line(OutPoints);
+	    dist = sqrt(dx * dx + dy * dy);
+	    Vect_append_point(OutPoints, Points->x[i - 1], Points->y[i - 1], Points->z[i - 1]);
+	}
+	Vect_append_point(OutPoints, Points->x[i], Points->y[i], Points->z[i]);
+    }
+    Vect_write_line(Map, otype, OutPoints, Cats);
+
+    Vect_destroy_line_struct(OutPoints);
+    /* G_free(OutPoints); */
+
+    return 0;
+}

vector/v.in.ogr/global.h

@@ -25,6 +25,7 @@
 
 
 extern int n_polygons;
+extern double split_distance;
 
 /* centroid structure */
 typedef struct

vector/v.in.ogr/main.c

@@ -38,6 +38,7 @@
 #endif
 
 int n_polygons;
+double split_distance;
 
 int geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
 	 double min_area, int type, int mk_centr);
@@ -55,9 +56,9 @@ int main(int argc, char *argv[])
 	*min_area_opt;
     struct Option *snap_opt, *type_opt, *outloc_opt, *cnames_opt;
     struct Flag *list_flag, *no_clean_flag, *z_flag, *notab_flag,
-	*region_flag;
+	*region_flag, *split_flag;
     struct Flag *over_flag, *extend_flag, *formats_flag, *tolower_flag;
-    char buf[2000], namebuf[2000];
+    char buf[2000], namebuf[2000], tempvect[GNAME_MAX];
     char *separator;
     struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
     struct Key_Value *proj_info, *proj_units;
@@ -65,7 +66,7 @@ int main(int argc, char *argv[])
     char error_msg[8192];
 
     /* Vector */
-    struct Map_info Map;
+    struct Map_info Map, Out;
     int cat;
 
     /* Attributes */
@@ -95,7 +96,13 @@ int main(int argc, char *argv[])
     int layer_id;
     int overwrite;
 
+    int feature, nfeatures;
+
+    double area_size = 0.;
+
     G_gisinit(argv[0]);
+ 
+    split_distance = -1.;
 
     module = G_define_module();
     module->keywords = _("vector, import");
@@ -233,6 +240,11 @@ int main(int argc, char *argv[])
 	_("Change column names to lowercase characters");
     tolower_flag->guisection = _("Attributes");
 
+    split_flag = G_define_flag();
+    split_flag->key = 's';
+    split_flag->description =
+	_("Split long boundaries to speed up cleaning");
+
     /* The parser checks if the map already exists in current mapset, this is
      * wrong if location options is used, so we switch out the check and do it
      * in the module after the parser */
@@ -452,6 +464,11 @@ int main(int argc, char *argv[])
 	cellhd.tb_res = 1.;
     }
 
+    if (split_flag->answer) {
+	split_distance = 0.;
+	area_size = sqrt((cellhd.east - cellhd.west) * (cellhd.north - cellhd.south));
+    }
+
     /* Fetch input map projection in GRASS form. */
     proj_info = NULL;
     proj_units = NULL;
@@ -581,14 +598,24 @@ int main(int argc, char *argv[])
     db_init_string(&sql);
     db_init_string(&strval);
 
-    /* open output vector */
+    /* open temporary vector, do the work in the temporary vector
+     * at the end copy alive lines to output vector
+     * in case of polygons this reduces the coor file size by a factor of 2 to 5
+     * only needed for polygons, but the presence of polygons can be detected
+     * only during OGR feature import, not before */
+    sprintf(buf, "%s", out_opt->answer);
+    /* strip any @mapset from vector output name */
+    G_find_vector(buf, G_mapset());
+    sprintf(tempvect, "%s_tmp", buf);
+    G_message(_("Using temporary vector <%s>"), tempvect);
     if (z_flag->answer)
-	Vect_open_new(&Map, out_opt->answer, 1);
+	Vect_open_new(&Map, tempvect, 1);
     else
-	Vect_open_new(&Map, out_opt->answer, 0);
+	Vect_open_new(&Map, tempvect, 0);
 
     Vect_hist_command(&Map);
 
+
     /* Points and lines are written immediately with categories. Boundaries of polygons are
      * written to the vector then cleaned and centroids are calculated for all areas in cleaan vector.
      * Then second pass through finds all centroids in each polygon feature and adds its category
@@ -778,9 +805,18 @@ int main(int argc, char *argv[])
 	cat = 1;
 	nogeom = 0;
 	OGR_L_ResetReading(Ogr_layer);
+	nfeatures = OGR_L_GetFeatureCount(Ogr_layer, 1);
+	feature = 0;
+	if (split_distance > -0.5 && nfeatures > 500) {
+	    split_distance = area_size / log(OGR_L_GetFeatureCount(Ogr_layer, 1));
+	    split_distance = split_distance / 3.; /* divisor is the handle */
+	    G_debug(1,"root of area size: %f", area_size);
+	    G_message(_("Boundary splitting distance in map units: %G"), split_distance);
+	}
 	G_important_message(_("Importing map %d features..."),
-			    OGR_L_GetFeatureCount(Ogr_layer, 1));
+			    nfeatures);
 	while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
+	    G_percent(feature++, nfeatures, 1);
 	    /* Geometry */
 	    Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
 	    if (Ogr_geometry == NULL) {
@@ -867,6 +903,7 @@ int main(int argc, char *argv[])
 	    OGR_F_Destroy(Ogr_feature);
 	    cat++;
 	}
+	G_percent(nfeatures, nfeatures, 1); /* finish it */
 
 	if (!notab_flag->answer) {
 	    db_commit_transaction(driver);
@@ -883,7 +920,7 @@ int main(int argc, char *argv[])
     G_message("%s", separator);
 
     /* TODO: is it necessary to build here? probably not, consumes time */
-    Vect_build(&Map);
+    Vect_build_partial(&Map, GV_BUILD_BASE);
 
     if (!no_clean_flag->answer &&
 	Vect_get_num_primitives(&Map, GV_BOUNDARY) > 0) {
@@ -898,12 +935,19 @@ int main(int argc, char *argv[])
 	Points = Vect_new_line_struct();
 
 	G_message("%s", separator);
-	G_warning(_("Cleaning polygons, result is not guaranteed!"));
 
+	/* first check if coor file size exceeds current 2GB limit */
+	/* coor file won't be opened if too large */
+	/* drawback: vector <tempvect> will persist as corrupted vector */
 	Vect_set_release_support(&Map);
 	Vect_close(&Map);
-	Vect_open_update(&Map, out_opt->answer, G_mapset());
-	Vect_build_partial(&Map, GV_BUILD_BASE);	/* Downgrade topo */
+	G_warning(_("Checking file size limits. On ERROR please delete vector <%s>"), tempvect);
+	Vect_open_update(&Map, tempvect, G_mapset());
+	/* topo was not fully built, build BASE again, is fast */
+	Vect_build_partial(&Map, GV_BUILD_BASE);
+
+	G_message("%s", separator);
+	G_warning(_("Cleaning polygons, result is not guaranteed!"));
 
 	if (snap >= 0) {
 	    G_message("%s", separator);
@@ -995,6 +1039,9 @@ int main(int argc, char *argv[])
 	    Vect_spatial_index_add_item(&si, centr, &box);
 	}
 
+	G_message("%s", separator);
+	G_message(_("Write centroids:"));
+
 	/* Go through all layers and find centroids for each polygon */
 	for (layer = 0; layer < nlayers; layer++) {
 	    G_message(_("Layer: %s"), layer_names[layer]);
@@ -1022,6 +1069,8 @@ int main(int argc, char *argv[])
 	for (centr = 1; centr <= ncentr; centr++) {
 	    double area;
 
+	    G_percent(centr, ncentr, 1);
+
 	    area = Vect_get_area_area(&Map, centr);
 	    total_area += area;
 
@@ -1052,11 +1101,11 @@ int main(int argc, char *argv[])
 	}
 	if (Centr)
 	    G_free(Centr);
-	G_message("%s", separator);
-	Vect_build_partial(&Map, GV_BUILD_NONE);
+	/* G_message("%s", separator); */
+	/* Vect_build_partial(&Map, GV_BUILD_NONE);
 
 	G_message("%s", separator);
-	Vect_build(&Map);
+	Vect_build(&Map); */
 
 	G_message("%s", separator);
 
@@ -1067,22 +1116,26 @@ int main(int argc, char *argv[])
 		      n_overlaps, nlayers + 1);
 	}
 
-	sprintf(buf, _("%d input polygons"), n_polygons);
-	G_message(buf);
+	Vect_hist_write(&Map, separator);
+	Vect_hist_write(&Map, "\n");
+	sprintf(buf, _("%d input polygons\n"), n_polygons);
+	G_message(_("%d input polygons"), n_polygons);
 	Vect_hist_write(&Map, buf);
 
-	sprintf(buf, _("Total area: %e (%d areas)"), total_area, ncentr);
-	G_message(buf);
+	sprintf(buf, _("Total area: %G (%d areas)\n"), total_area, ncentr);
+	G_message(_("Total area: %G (%d areas)"), total_area, ncentr);
 	Vect_hist_write(&Map, buf);
 
-	sprintf(buf, _("Overlapping area: %e (%d areas)"), overlap_area,
+	sprintf(buf, _("Overlapping area: %G (%d areas)\n"), overlap_area,
+		n_overlaps);
+	G_message(_("Overlapping area: %G (%d areas)"), overlap_area,
 		n_overlaps);
-	G_message(buf);
 	Vect_hist_write(&Map, buf);
 
-	sprintf(buf, _("Area without category: %e (%d areas)"), nocat_area,
+	sprintf(buf, _("Area without category: %G (%d areas)\n"), nocat_area,
+		n_nocat);
+	G_message(_("Area without category: %G (%d areas)"), nocat_area,
 		n_nocat);
-	G_message(buf);
 	Vect_hist_write(&Map, buf);
     }
 
@@ -1090,8 +1143,40 @@ int main(int argc, char *argv[])
      * OGR_DS_Destroy( Ogr_ds );
      */
 
+    /* second check if coor file size exceeds current 2GB limit */
+    /* coor file won't be opened if too large */
+    /* drawback: vector <tempvect> will persist as corrupted vector */
+    Vect_set_release_support(&Map);
     Vect_close(&Map);
+    G_message("%s", separator);
+    G_warning(_("Checking file size limits. On ERROR please delete vector <%s>"), tempvect);
+    Vect_open_old(&Map, tempvect, G_mapset());
 
+    /* Open output vector */
+    
+    /* Copy temporary vector to output vector */
+    G_message("%s", separator);
+    G_message(_("Copy temporary vector to output vector:"));
+    G_message("%s", separator);
+
+    if (z_flag->answer)
+	Vect_open_new(&Out, out_opt->answer, 1);
+    else
+	Vect_open_new(&Out, out_opt->answer, 0);
+
+    Vect_copy_head_data(&Map, &Out);
+    /* reduce coor file size by factor 2 to 5 in case of shapefile area import */
+    Vect_copy_map_lines(&Map, &Out);
+    Vect_copy_tables(&Map, &Out, 0);
+    Vect_hist_copy(&Map, &Out);
+
+    Vect_build(&Out);
+    
+    Vect_set_release_support(&Map);
+    Vect_set_release_support(&Out);
+    Vect_close(&Map);
+    Vect_close(&Out);
+    Vect_delete(tempvect);
 
     /* -------------------------------------------------------------------- */
     /*      Extend current window based on dataset.                         */

vector/v.net.iso/main.c

@@ -252,12 +252,13 @@ int main(int argc, char **argv)
 		continue;
 	    }			/* node unreachable */
 
+	    fprintf(stdout, "%f\n", cost);
+
 	    /* We must add centre node costs (not calculated by Vect_net_shortest_path() ), but
 	     *  only if centre and node are not identical, because at the end node cost is add later */
 	    if (node1 != node2)
 		cost += n1cost;
-	    G_debug(5,
-		    "Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
+	    G_debug(5, "Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
 		    node1, node2, cost, Nodes[node2].centre,
 		    Nodes[node2].cost);
 	    if (Nodes[node2].centre == -1 || cost < Nodes[node2].cost) {

vector/v.net.salesman/main.c

@@ -88,6 +88,7 @@ int main(int argc, char **argv)
     struct cat_list *Clist;
     struct line_cats *Cats;
     struct line_pnts *Points;
+    FILE *fp;
 
     /* Initialize the GIS calls */
     G_gisinit(argv[0]);
@@ -209,6 +210,8 @@ int main(int argc, char **argv)
     Vect_net_build_graph(&Map, type, afield, 0, afcol->answer, NULL, NULL,
 			 geo, 0);
 
+    fp = fopen("salesman_costs", "w+");
+
     /* Create sorted lists of costs */
     for (i = 0; i < ncities; i++) {
 	k = 0;
@@ -222,17 +225,22 @@ int main(int argc, char **argv)
 		G_fatal_error(_("Destination node [%d] is unreachable "
 				"from node [%d]"), cities[i], cities[j]);
 
+	    fprintf(fp, "%f\n", cost);
+
 	    costs[i][k].city = j;
 	    costs[i][k].cost = cost;
 	    k++;
 	}
 	qsort((void *)costs[i], k, sizeof(COST), cmp);
     }
+
+    fclose(fp);
+
     /* debug: print sorted */
     for (i = 0; i < ncities; i++) {
 	for (j = 0; j < ncities - 1; j++) {
 	    city = costs[i][j].city;
-	    G_debug(2, "%d -> %d = %f\n", cities[i], cities[city],
+	    G_debug(3, "%d -> %d = %f\n", cities[i], cities[city],
 		    costs[i][j].cost);
 	}
     }

vector/vectorintro.html

@@ -159,11 +159,11 @@ printed or maintained.
 It is possible to link the geographic objects in
 a vector map to one or more tables. Each link to a distinct
 attribute table is called a layer. A link defines which database
-driver, database and table is to be used. Each category number
-in a geometry file corresponds to a row in the attribute table
-(the linking column is usually the "cat" key column). Using
-<a href="v.db.connect.html">v.db.connect</a> layers can be listed
-or maintained.<br>
+driver, database and table is to be used. Each category number in a
+geometry file is associated with a layer and corresponds to a row in the
+attribute table for this layer (the linking column is usually the "cat"
+key column). Using <a href="v.db.connect.html">v.db.connect</a> layers
+can be listed or maintained.<br>
 Vector objects are not organized in layers. All vector
 objects are kept in one geometry file, and topology is maintained for
 all vector objects together. GRASS layers only consist of links to