[GRASS-SVN] r43710 -
grass/branches/releasebranch_6_4/vector/v.in.ogr
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Sep 27 14:03:56 EDT 2010
Author: mmetz
Date: 2010-09-27 18:03:56 +0000 (Mon, 27 Sep 2010)
New Revision: 43710
Modified:
grass/branches/releasebranch_6_4/vector/v.in.ogr/geom.c
grass/branches/releasebranch_6_4/vector/v.in.ogr/global.h
grass/branches/releasebranch_6_4/vector/v.in.ogr/main.c
Log:
Radim's TODO: split lines before breaking lines
Modified: grass/branches/releasebranch_6_4/vector/v.in.ogr/geom.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.in.ogr/geom.c 2010-09-27 18:02:01 UTC (rev 43709)
+++ grass/branches/releasebranch_6_4/vector/v.in.ogr/geom.c 2010-09-27 18:03:56 UTC (rev 43710)
@@ -25,6 +25,8 @@
#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
@@ -166,6 +168,35 @@
return 0;
}
+/* count polygons and isles */
+int poly_count(OGRGeometryH hGeom)
+{
+ int i, nr, ret;
+ OGRwkbGeometryType eType;
+ OGRGeometryH hRing;
+ eType = wkbFlatten(OGR_G_GetGeometryType(hGeom));
+
+ if (eType == wkbPolygon) {
+ G_debug(3, "Polygon");
+ nr = OGR_G_GetGeometryCount(hGeom);
+ n_polygon_boundaries += nr;
+
+ }
+ else if (eType == wkbGeometryCollection || eType == wkbMultiPolygon) {
+ G_debug(3, "GeometryCollection or MultiPolygon");
+ nr = OGR_G_GetGeometryCount(hGeom);
+ for (i = 0; i < nr; i++) {
+ hRing = OGR_G_GetGeometryRef(hGeom, i);
+
+ ret = poly_count(hRing);
+ if (ret == -1) {
+ G_warning(_("Cannot read part of geometry"));
+ }
+ }
+ }
+ return 0;
+}
+
/* Write geometry to output map */
int
geom(OGRGeometryH hGeom, struct Map_info *Map, int field, int cat,
@@ -223,7 +254,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) {
@@ -256,7 +291,7 @@
size = G_area_of_polygon(Points->x, Points->y, Points->n_points);
if (size < min_area) {
- G_warning(_("Area size [%.1e], area not imported"), size);
+ G_debug(2, "Area size [%.1e], area not imported", size);
return 0;
}
@@ -264,8 +299,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, BCats);
+ else
+ Vect_write_line(Map, otype, Points, BCats);
+
/* Isles */
IPoints =
(struct line_pnts **)G_malloc((nr - 1) *
@@ -298,7 +337,7 @@
IPoints[valid_isles]->y,
IPoints[valid_isles]->n_points);
if (size < min_area) {
- G_warning(_("Island size [%.1e], island not imported"),
+ G_debug(2, "Island size [%.1e], island not imported",
size);
}
else {
@@ -306,7 +345,11 @@
otype = GV_LINE;
else
otype = GV_BOUNDARY;
- Vect_write_line(Map, otype, IPoints[valid_isles], BCats);
+ if (split_distance > 0 && otype == GV_BOUNDARY)
+ split_line(Map, otype, IPoints[valid_isles], BCats);
+ else
+ Vect_write_line(Map, otype, IPoints[valid_isles],
+ BCats);
}
valid_isles++;
}
@@ -383,3 +426,52 @@
return 0;
}
+
+int split_line(struct Map_info *Map, int otype, struct line_pnts *Points,
+ struct line_cats *Cats)
+{
+ int i;
+ double dist = 0., seg_dist, dx, dy;
+ struct line_pnts *OutPoints;
+
+ /* don't write zero length boundaries */
+ Vect_line_prune(Points);
+ if (Points->n_points < 2)
+ return 0;
+
+ /* can't split boundaries with only 2 vertices */
+ if (Points->n_points == 2) {
+ Vect_write_line(Map, otype, Points, Cats);
+ return 0;
+ }
+
+ OutPoints = Vect_new_line_struct();
+ Vect_append_point(OutPoints, Points->x[0], Points->y[0], Points->z[0]);
+ Vect_append_point(OutPoints, Points->x[1], Points->y[1], Points->z[1]);
+ dx = Points->x[1] - Points->x[0];
+ dy = Points->y[1] - Points->y[0];
+ dist = sqrt(dx * dx + dy * dy);
+
+ /* trying to keep line length smaller than split_distance
+ * alternative, less code: write line as soon as split_distance is exceeded */
+ for (i = 2; i < Points->n_points; i++) {
+ dx = Points->x[i] - Points->x[i - 1];
+ dy = Points->y[i] - Points->y[i - 1];
+ seg_dist = sqrt(dx * dx + dy * dy);
+ dist += seg_dist;
+ if (dist > split_distance) {
+ Vect_write_line(Map, otype, OutPoints, Cats);
+ Vect_reset_line(OutPoints);
+ dist = seg_dist;
+ 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);
+
+ return 0;
+}
Modified: grass/branches/releasebranch_6_4/vector/v.in.ogr/global.h
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.in.ogr/global.h 2010-09-27 18:02:01 UTC (rev 43709)
+++ grass/branches/releasebranch_6_4/vector/v.in.ogr/global.h 2010-09-27 18:03:56 UTC (rev 43710)
@@ -26,8 +26,12 @@
#ifdef MAIN
int n_polygons = 0;
+int n_polygon_boundaries;
+double split_distance;
#else
extern int n_polygons;
+extern int n_polygon_boundaries;
+extern double split_distance;
#endif
Modified: grass/branches/releasebranch_6_4/vector/v.in.ogr/main.c
===================================================================
--- grass/branches/releasebranch_6_4/vector/v.in.ogr/main.c 2010-09-27 18:02:01 UTC (rev 43709)
+++ grass/branches/releasebranch_6_4/vector/v.in.ogr/main.c 2010-09-27 18:03:56 UTC (rev 43710)
@@ -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;
@@ -902,15 +966,10 @@
G_message("%s", separator);
G_warning(_("Cleaning polygons, result is not guaranteed!"));
- 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 */
-
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 +983,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 +1048,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,6 +1065,7 @@
/* Go through all layers and find centroids for each polygon */
for (layer = 0; layer < nlayers; layer++) {
+ G_message("%s", separator);
G_message(_("Layer: %s"), layer_names[layer]);
layer_id = layers[layer];
Ogr_layer = OGR_DS_GetLayer(Ogr_ds, layer_id);
@@ -1019,12 +1086,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 +1122,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 +1143,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