[GRASS-SVN] r43535 - grass/branches/develbranch_6/vector/v.in.ogr
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 20 04:13:29 EDT 2010
Author: mmetz
Date: 2010-09-20 08:13:28 +0000 (Mon, 20 Sep 2010)
New Revision: 43535
Modified:
grass/branches/develbranch_6/vector/v.in.ogr/main.c
Log:
activate boundary splitting
Modified: grass/branches/develbranch_6/vector/v.in.ogr/main.c
===================================================================
--- grass/branches/develbranch_6/vector/v.in.ogr/main.c 2010-09-20 08:04:55 UTC (rev 43534)
+++ grass/branches/develbranch_6/vector/v.in.ogr/main.c 2010-09-20 08:13:28 UTC (rev 43535)
@@ -41,6 +41,7 @@
double min_area, int type, int mk_centr);
int centroid(OGRGeometryH hGeom, CENTR * Centr, SPATIAL_INDEX * Sindex,
int field, int cat, double min_area, int type);
+int poly_count(OGRGeometryH hGeom);
int main(int argc, char *argv[])
{
@@ -55,7 +56,7 @@
struct Flag *list_flag, *no_clean_flag, *z_flag, *notab_flag,
*region_flag;
struct Flag *over_flag, *extend_flag, *formats_flag, *tolower_flag;
- char buf[2000], namebuf[2000];
+ char buf[2000], namebuf[2000], tempvect[2000];
char *separator;
struct Key_Value *loc_proj_info = NULL, *loc_proj_units = NULL;
struct Key_Value *proj_info, *proj_units;
@@ -63,7 +64,7 @@
char error_msg[8192];
/* Vector */
- struct Map_info Map;
+ struct Map_info Map, Tmp;
int cat;
/* Attributes */
@@ -92,6 +93,7 @@
int navailable_layers;
int layer_id;
int overwrite;
+ double area_size = 0.;
G_gisinit(argv[0]);
@@ -456,6 +458,16 @@
cellhd.tb_res = 1.;
}
+ /* split boundaries only when cleaning */
+ if (no_clean_flag->answer) {
+ split_distance = -1.;
+ }
+ else {
+ 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;
@@ -586,10 +598,24 @@
db_init_string(&strval);
/* open output vector */
- if (z_flag->answer)
+ /* 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_verbose_message(_("Using temporary vector <%s>"), tempvect);
+ if (z_flag->answer) {
Vect_open_new(&Map, out_opt->answer, 1);
- else
+ Vect_open_new(&Tmp, tempvect, 1);
+ }
+ else {
Vect_open_new(&Map, out_opt->answer, 0);
+ Vect_open_new(&Tmp, tempvect, 0);
+ }
Vect_hist_command(&Map);
@@ -622,7 +648,7 @@
if (ncnames > 0) {
cat_col_name = cnames_opt->answers[0];
}
- Vect_map_add_dblink(&Map, layer + 1, NULL, Fi->table,
+ Vect_map_add_dblink(&Map, layer + 1, layer_names[layer], Fi->table,
cat_col_name, Fi->database, Fi->driver);
ncols = OGR_FD_GetFieldCount(Ogr_featuredefn);
@@ -782,9 +808,45 @@
cat = 1;
nogeom = 0;
OGR_L_ResetReading(Ogr_layer);
- G_important_message(_("Importing map %d features..."),
- OGR_L_GetFeatureCount(Ogr_layer, 1));
+ unsigned int n_features = 0, feature_count = 0;
+
+ n_polygon_boundaries = 0;
+ n_features = OGR_L_GetFeatureCount(Ogr_layer, 1);
+
+ /* estimate distance for boundary splitting --> */
+
+ if (split_distance > -0.5) {
+ /* count polygons and isles */
+ G_message(_("Counting polygons for %d features..."), n_features);
+ while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
+ G_percent(feature_count++, n_features, 1); /* show something happens */
+ /* Geometry */
+ Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
+ if (Ogr_geometry != NULL) {
+ poly_count(Ogr_geometry);
+ }
+ OGR_F_Destroy(Ogr_feature);
+ }
+ /* rewind layer */
+ OGR_L_ResetReading(Ogr_layer);
+ feature_count = 0;
+ }
+
+ G_debug(1, "n polygon boundaries: %d", n_polygon_boundaries);
+ if (split_distance > -0.5 && n_polygon_boundaries > 50) {
+ split_distance =
+ area_size / log(n_polygon_boundaries);
+ /* divisor is the handle: increase divisor to decrease split_distance */
+ split_distance = split_distance / 4.;
+ G_debug(1, "root of area size: %f", area_size);
+ G_verbose_message(_("Boundary splitting distance in map units: %G"),
+ split_distance);
+ }
+ /* <-- estimate distance for boundary splitting */
+
+ G_important_message(_("Importing map %d features..."), n_features);
while ((Ogr_feature = OGR_L_GetNextFeature(Ogr_layer)) != NULL) {
+ G_percent(feature_count++, n_features, 1); /* show something happens */
/* Geometry */
Ogr_geometry = OGR_F_GetGeometryRef(Ogr_feature);
if (Ogr_geometry == NULL) {
@@ -795,7 +857,7 @@
if (dim > 2)
with_z = 1;
- geom(Ogr_geometry, &Map, layer + 1, cat, min_area, type,
+ geom(Ogr_geometry, &Tmp, layer + 1, cat, min_area, type,
no_clean_flag->answer);
}
@@ -869,6 +931,7 @@
OGR_F_Destroy(Ogr_feature);
cat++;
}
+ G_percent(n_features, n_features, 1); /* finish it */
if (!notab_flag->answer) {
db_commit_transaction(driver);
@@ -885,10 +948,11 @@
G_message("%s", separator);
/* TODO: is it necessary to build here? probably not, consumes time */
- Vect_build(&Map);
+ /* GV_BUILD_BASE is sufficient to toggle boundary cleaning */
+ Vect_build_partial(&Tmp, GV_BUILD_BASE);
if (!no_clean_flag->answer &&
- Vect_get_num_primitives(&Map, GV_BOUNDARY) > 0) {
+ Vect_get_num_primitives(&Tmp, GV_BOUNDARY) > 0) {
int ret, centr, ncentr, otype, n_overlaps, n_nocat;
CENTR *Centr;
SPATIAL_INDEX si;
@@ -910,7 +974,7 @@
if (snap >= 0) {
G_message("%s", separator);
G_message(_("Snap boundaries (threshold = %.3e):"), snap);
- Vect_snap_lines(&Map, GV_BOUNDARY, snap, NULL);
+ Vect_snap_lines(&Tmp, GV_BOUNDARY, snap, NULL);
}
/* It is not to clean to snap centroids, but I have seen data with 2 duplicate polygons
@@ -924,57 +988,64 @@
G_message("%s", separator);
G_message(_("Break polygons:"));
- Vect_break_polygons(&Map, GV_BOUNDARY, NULL);
+ Vect_break_polygons(&Tmp, GV_BOUNDARY, NULL);
- /* It is important to remove also duplicate centroids in case of duplicate imput polygons */
+ /* It is important to remove also duplicate centroids in case of duplicate input polygons */
G_message("%s", separator);
G_message(_("Remove duplicates:"));
- Vect_remove_duplicates(&Map, GV_BOUNDARY | GV_CENTROID, NULL);
+ Vect_remove_duplicates(&Tmp, GV_BOUNDARY | GV_CENTROID, NULL);
+ /* split boundaries here ? */
+
/* Vect_clean_small_angles_at_nodes() can change the geometry so that new intersections
* are created. We must call Vect_break_lines(), Vect_remove_duplicates()
- * and Vect_clean_small_angles_at_nodes() until no more small dangles are found */
+ * and Vect_clean_small_angles_at_nodes() until no more small angles are found */
do {
G_message("%s", separator);
G_message(_("Break boundaries:"));
- Vect_break_lines(&Map, GV_BOUNDARY, NULL);
+ Vect_break_lines(&Tmp, GV_BOUNDARY, NULL);
G_message("%s", separator);
G_message(_("Remove duplicates:"));
- Vect_remove_duplicates(&Map, GV_BOUNDARY, NULL);
+ Vect_remove_duplicates(&Tmp, GV_BOUNDARY, NULL);
G_message("%s", separator);
G_message(_("Clean boundaries at nodes:"));
nmodif =
- Vect_clean_small_angles_at_nodes(&Map, GV_BOUNDARY, NULL);
+ Vect_clean_small_angles_at_nodes(&Tmp, GV_BOUNDARY, NULL);
} while (nmodif > 0);
G_message("%s", separator);
if (type & GV_BOUNDARY) { /* that means lines were converted boundaries */
G_message(_("Change boundary dangles to lines:"));
- Vect_chtype_dangles(&Map, -1.0, NULL);
+ Vect_chtype_dangles(&Tmp, -1.0, NULL);
}
else {
G_message(_("Change dangles to lines:"));
- Vect_remove_dangles(&Map, GV_BOUNDARY, -1.0, NULL);
+ Vect_remove_dangles(&Tmp, GV_BOUNDARY, -1.0, NULL);
}
G_message("%s", separator);
if (type & GV_BOUNDARY) {
G_message(_("Change boundary bridges to lines:"));
- Vect_chtype_bridges(&Map, NULL);
+ Vect_chtype_bridges(&Tmp, NULL);
}
else {
G_message(_("Remove bridges:"));
- Vect_remove_bridges(&Map, NULL);
+ Vect_remove_bridges(&Tmp, NULL);
}
+ /* merge boundaries */
+ G_message("%s", separator);
+ G_message(_("Merge boundaries:"));
+ Vect_merge_lines(&Tmp, GV_BOUNDARY, NULL, NULL);
+
/* Boundaries are hopefully clean, build areas */
G_message("%s", separator);
- Vect_build_partial(&Map, GV_BUILD_ATTACH_ISLES);
+ Vect_build_partial(&Tmp, GV_BUILD_ATTACH_ISLES);
/* Calculate new centroids for all areas, centroids have the same id as area */
- ncentr = Vect_get_num_areas(&Map);
+ ncentr = Vect_get_num_areas(&Tmp);
G_debug(3, "%d centroids/areas", ncentr);
Centr = (CENTR *) G_calloc(ncentr + 1, sizeof(CENTR));
@@ -982,7 +1053,7 @@
for (centr = 1; centr <= ncentr; centr++) {
Centr[centr].valid = 0;
Centr[centr].cats = Vect_new_cats_struct();
- ret = Vect_get_point_in_area(&Map, centr, &x, &y);
+ ret = Vect_get_point_in_area(&Tmp, centr, &x, &y);
if (ret < 0) {
G_warning(_("Cannot calculate area centroid"));
continue;
@@ -999,7 +1070,8 @@
/* Go through all layers and find centroids for each polygon */
for (layer = 0; layer < nlayers; layer++) {
- G_message(_("Layer: %s"), layer_names[layer]);
+ G_message("%s", separator);
+ G_message(_("Find centroids for layer: %s"), layer_names[layer]);
layer_id = layers[layer];
Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layer_id);
OGR_L_ResetReading(Ogr_layer);
@@ -1019,12 +1091,17 @@
}
/* Write centroids */
+ G_message("%s", separator);
+ G_message(_("Write centroids:"));
+
n_overlaps = n_nocat = 0;
total_area = overlap_area = nocat_area = 0.0;
for (centr = 1; centr <= ncentr; centr++) {
double area;
+
+ G_percent(centr, ncentr, 2);
- area = Vect_get_area_area(&Map, centr);
+ area = Vect_get_area_area(&Tmp, centr);
total_area += area;
if (!(Centr[centr].valid)) {
@@ -1050,17 +1127,19 @@
otype = GV_POINT;
else
otype = GV_CENTROID;
- Vect_write_line(&Map, otype, Points, Centr[centr].cats);
+ Vect_write_line(&Tmp, otype, Points, Centr[centr].cats);
}
if (Centr)
G_free(Centr);
- G_message("%s", separator);
- Vect_build_partial(&Map, GV_BUILD_NONE);
+
+ Vect_spatial_index_destroy(&si);
- G_message("%s", separator);
- Vect_build(&Map);
+ /* G_message("%s", separator); */
+ /* Vect_build_partial(&Map, GV_BUILD_NONE); */
- G_message("%s", separator);
+ /* not needed */
+ /* G_message("%s", separator);
+ Vect_build(&Map); */
if (n_overlaps > 0) {
G_warning(_("%d areas represent more (overlapping) features, because polygons overlap "
@@ -1069,29 +1148,42 @@
n_overlaps, nlayers + 1);
}
- sprintf(buf, _("%d input polygons"), n_polygons);
- G_message(buf);
+ G_message("%s", separator);
+
+ 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);
+ G_message("%s", separator);
}
/* needed?
* OGR_DS_Destroy( Ogr_ds );
*/
+ /* Copy temporary vector to output vector */
+ Vect_copy_map_lines(&Tmp, &Map);
+ Vect_close(&Tmp);
+ Vect_delete(tempvect);
+
+ Vect_build(&Map);
Vect_close(&Map);
More information about the grass-commit
mailing list