[GRASS-SVN] r36310 - in grass/trunk/vector: . v.in.ogr v.net.iso v.net.salesman

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Mar 10 10:00:32 EDT 2009


Author: mmetz
Date: 2009-03-10 10:00:32 -0400 (Tue, 10 Mar 2009)
New Revision: 36310

Modified:
   grass/trunk/vector/v.in.ogr/geom.c
   grass/trunk/vector/v.in.ogr/global.h
   grass/trunk/vector/v.in.ogr/main.c
   grass/trunk/vector/v.net.iso/main.c
   grass/trunk/vector/v.net.salesman/main.c
   grass/trunk/vector/vectorintro.html
Log:
vectorintro corrected

Modified: grass/trunk/vector/v.in.ogr/geom.c
===================================================================
--- grass/trunk/vector/v.in.ogr/geom.c	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/v.in.ogr/geom.c	2009-03-10 14:00:32 UTC (rev 36310)
@@ -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 @@
 	    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,8 +248,12 @@
 	    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 =
 	    (struct line_pnts **)G_malloc((nr - 1) *
@@ -274,7 +283,11 @@
 		    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 @@
 
     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;
+}

Modified: grass/trunk/vector/v.in.ogr/global.h
===================================================================
--- grass/trunk/vector/v.in.ogr/global.h	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/v.in.ogr/global.h	2009-03-10 14:00:32 UTC (rev 36310)
@@ -25,6 +25,7 @@
 
 
 extern int n_polygons;
+extern double split_distance;
 
 /* centroid structure */
 typedef struct

Modified: grass/trunk/vector/v.in.ogr/main.c
===================================================================
--- grass/trunk/vector/v.in.ogr/main.c	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/v.in.ogr/main.c	2009-03-10 14:00:32 UTC (rev 36310)
@@ -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 @@
 	*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 @@
     char error_msg[8192];
 
     /* Vector */
-    struct Map_info Map;
+    struct Map_info Map, Out;
     int cat;
 
     /* Attributes */
@@ -95,7 +96,13 @@
     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 @@
 	_("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 @@
 	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 @@
     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 @@
 	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 @@
 	    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 @@
     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,13 +935,20 @@
 	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);
 	    G_message(_("Snap boundaries (threshold = %.3e):"), snap);
@@ -995,6 +1039,9 @@
 	    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 @@
 	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 @@
 	}
 	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 @@
 		      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(buf);
+	G_message(_("Overlapping area: %G (%d areas)"), overlap_area,
+		n_overlaps);
 	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(buf);
+	G_message(_("Area without category: %G (%d areas)"), nocat_area,
+		n_nocat);
 	Vect_hist_write(&Map, buf);
     }
 
@@ -1090,9 +1143,41 @@
      * 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.                         */
     /* -------------------------------------------------------------------- */

Modified: grass/trunk/vector/v.net.iso/main.c
===================================================================
--- grass/trunk/vector/v.net.iso/main.c	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/v.net.iso/main.c	2009-03-10 14:00:32 UTC (rev 36310)
@@ -252,12 +252,13 @@
 		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) {

Modified: grass/trunk/vector/v.net.salesman/main.c
===================================================================
--- grass/trunk/vector/v.net.salesman/main.c	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/v.net.salesman/main.c	2009-03-10 14:00:32 UTC (rev 36310)
@@ -88,6 +88,7 @@
     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 @@
     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 @@
 		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);
 	}
     }

Modified: grass/trunk/vector/vectorintro.html
===================================================================
--- grass/trunk/vector/vectorintro.html	2009-03-10 13:21:19 UTC (rev 36309)
+++ grass/trunk/vector/vectorintro.html	2009-03-10 14:00:32 UTC (rev 36310)
@@ -159,11 +159,11 @@
 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



More information about the grass-commit mailing list